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ABSTRACT 

Providing  a  simplified  model  of  real  terrain  has  applications  to  route  planning  for 
robotic  vehicles  and  militarv'  maneuvers.  In  this  thesis  I  explore  planar-patch  surface 
modeling  to  represent  terrain  in  a  simple  and  effective  way.  In  planar-patch  surface 
modeling  the  terrain  is  subdivided  into  a  set  of  planar  subregions.  The  homogeneity  of 
the  gradient  within  a  planar  subregion  simplifies  calculating  the  cost  of  traversing  the 
region,  thus  simplifying  route  planning.  I  have  explored  three  main  strategies  to  model 
the  surface:  joint  top-down  and  and  bottom-up,  strict  bottom-up,  and  presmootliing 
bottom-up  approaches.  Results  of  the  algorithms  are  shown  graphically  by  using  the 
APL  and  Grafstat  packages,  verifying  their  correctness  and  accuracy. 


m 


TABLE  OF  CONTENTS 

I.  INTRODUCTION   1 

A.  BACKGROUND    1 

B.  ORGANIZATION    1 

II.  REVIEW  OF  TERRAIN  MODELING     3 

A.  TERRAIN  CELL  CLASSIFICATION    3 

1.  QUADR^^TIC  APPROXIMATION  OF  SURFACE  BY  LEAST 
SQUARES  FIT    3 

2.  TERRAIN  CELL  CLASSIFICATION  BASED  ON  QUADRATIC  AP- 
PROXIMATIONS     5 

B.  PLANE  APPROXIMATION  OF  A  SURFACE  BY  LEAST  SQUARES  FIT   6 

C.  THE  GRADIENT  AT  A  POINT     7 

III.  PLANAR-PATCH  TERRAIN  MODELS     11 

A.  IMPLEMENTATION     II 

B.  JOINT  TOP-DOWN  AND  BOTTOM-UP  PLANAR-PATCH  TERRAIN 
MODELING    12 

1.  TOP-DOWN  PHASE    13 

2.  THRESHOLDS  USED   13 

3.  BOTTOM-UP  PHASE    15 

C.  STRICT  BOTTOM-UP  PLANAR-PATCH  TERRAIN  MODELING    18 

1.  POINT-CLASSIFICATION  PHASE     18 

2.  REGION-GROWING  PHASE   18 

a.  ONE-DIMENSIONAL  REGION-GROWING     18 

b.  TWO-DIMENSIONAL  REGION-GROWING    19 

D.  AN  ALGORITHM  IMPROVING  THE  CONTINUITY  OF  ADJOINING 
REGIONS    21 

IV.  DISCUSSION    24 

A.  ACCURACY  ANALYSIS     24 

B.  LIMITATIONS  OF  THE  MODEL   24 


IV 


■  LIEB  A  RY 

C.     EXPERIMENTAL  RESULTS     25 

L     EXPERIMENTAL  RESULTS  1  OF  THE  FIRST  APPROACH     26 

2.  EXPERIMENTAL  RESULTS  2  OF  THE  FIRST  APPROACH     31 

3.  EXPERIMENTAL  RESULTS  1  OF  THE  SECOND  APPROACH    ...  35 

4.  EXPERIMENTAL  RESULTS  2  OF  THE  SECOND  APPROACH    ...  39 

5.  EXPERIMENTAL  RESULTS  OF  THE  THIRD  APPROACH    43 

V.     CONCLUSIONS   ' 48 

APPENDIX  A.     PASCAL  CODE  OF  JOINT  TOP-DOWN  AND  BOTTOM-UP 
ALGORITHM    49 

APPENDIX  B.     PASCAL  CODE  OF  STRICT  BOTTOM-UP  ALGORITHM    ...  66 

APPENDIX  C.     PASCAL  CODE  FOR  SMOOTHING  ELEVATION  DATA  AND 
CALCULATING  GR^^DIENT     85 

LIST  OF  REFERENCES   90 

INITIAL  DISTRIBUTION  LIST     91 


LIST  OF  TABLES 

Table     1.  CLASSIFICATION  OF  SURFACE  POINTS  BY  EIGENVALUFS  OF 

THE  HESSIAN  MATRIX    6 

Table     2.  A  MATRIX  OF  TERRAIN  POINTS     8 

Table     3.  AN  EXAMPLE  OF  A  TERR.-\IN  POINTS  MATRIX    9 


VI 


LIST  OF  FIGURES 

Figure  1.  Block  Diagram  of  the  First  Two  Program  Approaches    12 

Figure  2.  Quadtree  Region  Subdivision  Method    14 

Figure  3.  Pascal  Code  for  the  Quadtree  Region  Subdivision  Method    15 

Figure  4.  Minimum-Size  Subregions  Created  by  the  Quadtree  Method     16 

Figure  5.  Subregions  after  the  Merging  Process     17 

Figure  6.  One-Dimensional  Region-Growing     19 

Figure  7.  Two-Dimensional  Region-Growing    20 

Figure  8.  Filling  the  Boundary  with  Two  Triangles    22 

Figure  9.  How  Surface  Convexity  Helps    23 

Figure  10.  Sample  Real  Terrain    26 

Figure  11.  Subregions  Created  by  the  Quadtree  Method  in  the  First  Approach    ...  27 

Figure  12.  Subregions  Created  after  the  Merging  in  the  First  Approach     28 

Figure  13.  Terrain  Model  Created  by  the  First  Approach    29 

Figure  14.  Terrain  Model  Created  by  the  First  Approach  Showing  Boundaries     ...  30 

Figure  15.  Subregions  Created  by  the  Quadtree  Method  in  the  First  Approach    ...  31 

Figure  16.  Subregions  Created  after  the  Merging  in  the  First  Approach     32 

Figure  17.  Terrain  Model  Created  by  the  First  Approach    33 

Figure  18.  Terrain  Model  Created  by  the  First  Approach  Showing  Boundaries     ...  34 

Figure  19.  Subregions  Created  by  the  Region  Growing  Algorithm  in  the  Second 

Approach    35 

Figure  20.  Subregions  Created  after  the  Merging  in  the  Second  Approach     36 

Figure  21.  Terrain  Model  Created  by  the  Second  Approach    37 

Figure  22.  Terrain  Model  Created  by  the  Second  Approach  Showing  Boundaries     .  38 

Figure  23.  Subregions  Created  by  the  Region  Growing  Algorithm  in  the  Second 

Approach    39 

Figure  24.  Subregions  Created  after  the  Merging  in  the  Second  Approach     40 

Figure  25.  Terrain  Model  Created  by  the  Second  Approach    41 

Figure  26.  Terrain  Model  Created  by  the  Second  Approach  Showing  Boundaries     .  42 

Figure  27.  Smoothed  Sample  Terrain    43 

Figure  28.  Subregions  Created  by  the  Third  Approach    44 

Figure  29.  Subregions  Created  after  the  Merging  in  the  Third  Approach    45 


vu 


Figure  30.  Terrain  Model  Created  by  the  Third  Approach     46 

Figure  31.  Terrain  Model  Created  by  the  Third  Approach  Showing  Boundaries    ...  47 


vm 


I.     INTRODUCTION 

A.  BACKGROUND 

There  are  many  ways  to  represent  a  three-dimensional  object,  usually  involving 
complicated  mathematical  and  graphical  theories.  Representing  and  displaying  a  three- 
dimensional  object  like  a  terrain  is  important  in  applications  such  as  cartography,  route 
planning  for  robotic  vehicles,  road  construction,  and  planning  of  military  maneuvers. 
A  simpUfied  model  of  real  terrain  helps  the  user  better  understand  and  analyze  it. 

The  simplest  approximations  are  piecewise  planar,  and  the  very  simplest  is  a 
polyhedron  with  triangular  faces.  Assume  that  a  terrain  is  approximated  by  a  two- 
dimensional  array  of  evenly  spaced  elevations.  Each  group  of  three  adjacent  points  in 
the  array  defines  a  planar  polyhedral  face,  and  a  set  of  triplets  constitutes  a  planar-patch 
approximation  of  the  surface.  Unfortunately,  this  approach  creates  many  small  triangu- 
lar regions,  which  can  make  the  searching  process  of  route  planning  too  complicated. 

In  this  thesis,  I  developed  three  algorithms  for  planar-patch  modeling  of  evenly 
spaced  elevation  data,  and  each  has  its  own  advantages  and  disadvantages.  I  have  used 
the  least-squares  method  to  fit  a  plane  over  the  points  of  a  terrain,  a  quadtree  algorithm 
to  subdivide  regions,  a  gradient  clustering  method  to  build  regions,  and  smoothing  al- 
gorithms on  elevation  data  to  improve  continuity  of  adjoining  regions. 

One  major  advantage  of  planar-patch  modeling  over  other  three-dimensional  surface 
modeling  is  that  terrain  within  a  patch  is  homogeneous  in  the  sense  of  a  constant  gra- 
dient. This  homogeneity  simplifies  calculating  the  cost  or  speed  of  traversing  the  region, 
thus  making  it  easier  to  analyze  for  route  planning. 

Several  thresholds  affect  my  algorithms,  and  different  threshold  values  will  result  in 
difierent  models.  These  include  the  standard  deviation  of  a  fitted  plane  from  given 
points,  the  ratio  of  magnitude  between  adjacent  points,  and  the  root  mean  square  dif- 
ference of  coefiicients  of  planar  equations  between  adjacent  regions.  The  final  result  of 
this  thesis  will  be  provided  to  Ron  Ross  (Ph.D.  student.  Computer  Science  Department, 
Naval  Postgraduate  School)  for  his  research  in  route  planning  for  robotic  vehicles. 

B.  ORGANIZATION 

Chapter  2  introduces  terrain  modeling  techniques,  including  least-squares  planar 
fittin2  and  cell  classification. 


Chapter  3  is  the  main  chapter  of  this  thesis,  defining  the  planar-patch  terrain  model. 
Joint  top-down  and  bottom-up  and  strict  bottom-up  planar-patch  terrain  modeling 
strategies  are  described  in  detail.  Programming  techniques  such  as  the  quadtree  method 
used  in  the  joint  top-down  and  bottom-up  algorithm  and  the  region-growing  method 
used  in  the  strict  bottom-up  algorithm  are  explained.  Finally,  the  issue  of  continuity  of 
adjacent  planar-patches  along  their  common  boundaries  is  discussed.  Attempts  to  solve 
this  problem,  such  as  smoothing  of  original  elevation  data  and  grouping  of  the  same 
type  of  points  into  one  region  are  explained. 

In  Chapter  4  comparisons  are  made  between  the  two  algorithms  in  terms  of  accu- 
racy and  simplicity  of  terrain  representation.  Special  features  and  the  limitations  inher- 
ent to  each  algorithm  are  discussed.  Experimental  results  of  joint  top-down  and 
bottom-up,  strict  bottom-up,  and  presmoothing  bottom-up  algorithms  are  presented  in 
graphical  form  in  this  chapter.  Pictures  of  the  modeled  terrain  are  shown  in  three  di- 
mensions using  pictures  drawn  by  the  Grafstat  and  APL  packages.  Pictures  of  the 
boundaries  of  subregions  are  also  included. 

Finally,  Chapter  5  summarizes  the  contributions  made  by  this  research  and  discusses 
some  of  the  possible  areas  for  further  research  based  on  this  work.  The  Pascal  source 
program  is  in  the  Appendix. 


II.     REVIEW  OF  TERRAIN  MODELING 

The  most  common  method  of  terrain  representation  is  by  a  set  of  functional  ex- 
pressions such  as 

J{x^-)  =  x   +3xy-5y+y 

where  f(x,y)  is  the  elevation  function  for  terrain  and 
X  and  y  are  the  latitude  and  longitude  coordinates. 

This  permits  the  simple  calculation  of  the  elevation  given  any  combination  of  x  and  y. 
It  is  known  that  any  continuous  surface  can  be  approximated  with  an  arbitrarily  small 
error  by  a  polynomial  of  sufficiently  high  degree.  Such  polynomials  can  be  derived  by 
least-squares  surface-fitting  methods. 

Surfaces  may  be  defined  in  tabular  form  where  the  value  of  f  is  given  for  a  selection 
of  representative  arguments.  To  find  an  arbitrary  value  off,  a  table  lookup  can  be  per- 
formed to  determine  the  value  of  the  nearest  known  points  and  use  them  to  approximate 
the  desired  value  by  interpolation.   [Ref  1:  p.  212] 

It  is  possible  to  specify  a  surface  in  terms  of  guiding  points  by  a  generalization  of 
the  Bezier  polynomials  or  the  B-splines.  The  Bezier  polynomial  or  B-splinc  method  finds 
an  approximate  curve  that  passes  near  a  set  of  given  points.   [Ref  2:  p.  247] 

A.     TERRAIN  CELL  CLASSIFICATION 

I.     QUADRATIC  APPROXIMATION  OF  SURFACE  BY  LEAST  SQUARES 
FIT 

Fitting  a  least-squares  quadratic  to  a  3  by  3  square  of  points  on  a  surface  is  done  in  the 
following  way.   The  fit  as  a  function  of  x  and  y  is:  [Ref  3:  p.  55] 

A  2  2 

J{x^)  =  k^  +  k2X  +  k^y  +  k^x   +  k^xy  +  k^ 

The  square  of  the  error  between  the  fit  equation  and  the  given  terrain  points  is: 

E^  =  2_jZj(^\  +  ^2-^  +  ^3>'  +  k^x^k^xy  +  k^^  -Ax^)f 

X     y 

The  partial  derivatives  of  the  squared  error  with  respect  to  each  constant  are: 
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To  minimize  the  squared  error,  each  of  the  above  equations  is  set  to  zero  and  simplified. 
Those  equations  can  be  further  simplified  by  assuming  that  the  coordinates  of  each  of 
the  points  in  the  3  by  3  matrix,  except  the  center  point,  is  ±  one  distance  unit  from  the 
center  point.  By  the  above  assumption,  any  terms  in  above  equations  which  contain  a 
summation  of  an  x  or  y  coordinate  with  an  odd  power  will  go  to  zero.  Simplified 
equations  are: 


0  =  A-,y  y  +  A4  X  6  +  A^  X  6  -  ^^.xrvV) 


(11) 
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0  =  A,  x6  +  A4x4  +  A6^6-  ^^X/^-^JOj'^  (16) 

From  above  equations,  if  we  solve  for  Aj,  Aj,  A3,  A4,  A5  and  A^,  we  get  a  quadratic 
equation  fitted  to  the  3  by  3  square  of  points. 

2.     TERRAIN  CELL  CLASSIFICATION  BASED  ON  QUADRATIC 
APPROXIMATIONS 

The  eigenvalues  of  the  Hessian  matrix  for  the  quadratic  function  approximated 
for  the  3  by  3  points  on  the  surface  are:  [Ref  3:  p.  61] 


(2A,  +  2k,)  ±  v' (  -2A,  -  2A,)'  -  4(2A42A,  -  A^) 
/  = ;: 


Using  the  signs  of  the  eigenvalues,  we  can  classify  the  central  point  by  Table  1 
on  page  6.  [Ref  3:  p.  69] 


Table  1.  CLASSIFICATION  OF  SURFACE 
POINTS  BY  EIGENVALUES  OF  THE 
HESSIAN  MATRIX 


sign  of/. 

sign  of /^.2 

- 

0 

+ 

+ 

saddle 

impossible 

depression 

0 

ridge 

plane 

valley 

- 

hill 

impossible 

pass 

B.     PLANE  APPROXIMATION  OF  A  SURFACE  BY  LEAST  SQUARES  FIT 

The  elevation  z  of  any  point  (x,y)  can  be  expressed  as: 

The  plane  that  we  want  to  fit  over  the  surface  is: 

z  =  ax  +  by  +  c 

The  square  of  the  error  between  the  fitted  plane  and  the  real  terrain  is: 


E^  =  XZ^""^  +  ^>'  "^  ^  --A^vV-))^ 


X     y 


The  partial  derivatives  of  the  squared  error  w'ith  respect  to  a,b  and  c  are: 
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To  minimize  the  squared  error,  each  of  the  above  equations  is  expanded  and  then  the 
partials  are  set  to  zero. 
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Here,  we  have  three  unknowns,  a,b  and  c,  and  three  equations.  So  we  can  solve  for 
them. 

C.     THE  GRADIENT  AT  A  POINT 

One  of  the  primar}'  criterion  we  will  use  for  grouping  points  into  a  region  is  the 
gradient.  The  gradient  of  a  point  is  a  vector  that  has  a  magnitude  (slope)  and  a  direction 
for  its  components.  We  can  represent  a  three  by  three  elevation  matrix  of  terrain  points 
as  in  Table  2  on  page  8,  where  x  and  y  are  the  coordinates  of  a  point  on  a  terrain  and 
fl^x,y)  is  the  elevation  of  the  point. 


Table  2.     A   MATRIX   OF   TERRAIN 
POINTS 


fix-l.y+l) 

flx,y+l) 

flx+l,y+l) 

f(x-I,y) 

flx.y) 

flx+l,y) 

r(x-l,y-l) 

flx,y-l) 

Rx+l,y-l) 

The  gradient  vector  of  a  point  {x,y)  can  be  approximated  in  three  ways.    The  first 
method  is: 

A,  =J{x+  \^v)-f[x^) 


l(\^  A.  \h 


magni[ude{x^v)  =  ^  (A,  +  A2) 


-1    ^2 
direciionixAi)  =  tan    (  — —  ) 

A] 


The  second  method  is: 


A,  =/(^+ij-)-/(^-ivy) 


^2=Ax^^+\)-Axy-\) 


The  magnitude(x,y)  and  direction(x,y)  are  the  same  as  above. 


The  third  way  of  approximating  gradient  is  by  first  fitting  a  quadratic  to  the  three 
by  three  square  of  points  and  then  using  the  coefficients  of  the  quadratic  :  [Ref  3:  p.  60] 


magnitude{.x^-)  =  >^/  (A'f  +  k^) 


-1.^3 


direction{x^)  =  tan    (  — —  ) 


where  /c^  and  A'3  are  the  coefficients  of  the  quadratic  described  in  section  A. 

The  first  method  considers  only  one  quadrant  of  the  3  by  3  matrix,  thus  is  biased. 
The  second  method  uses  a  point  from  each  of  the  quadrants,  and  is  better  than  the  first 
method  in  the  sense  that  it  utiUzes  more  unbiased  information.  However,  it  gives  an 
incorrect  value  in  the  case  as  shown  in  Table  3. 


Table  3.     AN  EXAMPLE  OF  A  TER- 
RAIN POINTS  MATRIX 


5 

2 

3 

3 

1 

3 

3 

2 

4 

When  we  calculate  the  gradient  of  center  point  using  the  second  method,  it  will  give 
the  magnitude  value  zero  and  direction  value  45  degree,  which  is  wrong.  The  third 
method  does  not  have  such  problem,  since  it  utilizes  all  eight  neighboring  points.  In  this 
case  it  gives  the  magnitude  value  0.2357  and  direction  value  45  degree.  I  used  the  third 
method  in  my  program. 


The  gradient  vector  of  a  terrain  point  points  in  the  x-y  direction  of  maximum  slope 
of  that  point;  magnitude'^ x,y)  is  the  maximum  slope  value  and  direction(x,y)  is  the  di- 
rection of  that  slope  in  the  x-y  plane. 
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III.     PLANAR-PATCH  TERRAIN  MODELS 

A  planar-patch  terrain  model  is  a  model  that  represents  terrain  in  the  form  of  a 
polyhedron.  This  model  is  different  from  the  common  planar  terrain  model  used  in 
graphics  in  that  size  of  a  planar-patch  is  not  fixed,  and  the  shape  of  a  patch  does  not 
have  to  be  regular. 

A.     IMPLEMENTATION 

I  chose  Pascal  as  an  implementation  language  for  this  research  because  of  my  fa- 
miliarity and  experience,  and  the  Waterloo  Pascal  on  the  IBM  370/3033AP  mainframe 
was  the  Pascal  compiler  that  was  used. 

The  APL  language  and  the  Grafstat  graphics  software  package  were  used  to  verify 
the  correctness  and  accuracy  of  the  algorithms.  Even  though  I  chose  Pascal  as  my  im- 
plementation language,  the  matrix  form  of  the  elevation  data  suggests  APL.  With  the 
Grafstat  software  package  which  runs  in  the  APL  environment,  one  can  check  the 
intermediate  results  of  an  APL  program  by  displaying  pictures.  Since  I  experimented 
with  different  threshold  values  in  my  algorithm,  APL  and  Grafstat  were  useful.  The 
three-dimensional  pictures  in  chapter  4  were  done  with  Grafstat. 

Two  arrays  of  records  are  used  as  main  data  structures  by  all  my  algorithms: 
point_array  and  subregion_array.  The  point_array  is  a  two-dimensional  array  which 
corresponds  to  the  points  m  the  real  terrain.  Each  element  of  the  array  is  a  record  which 
has  the  elevation,  the  subregion  ID,  and  the  magnitude  of  the  terrain  point  as  its  fields. 
From  the  values  of  elevation  and  magnitude,  the  subregion  ID  is  determined  and  as- 
signed to  the  point.  The  subregion__array  is  a  one-dimensional  array  of  records  which  is 
created  during  the  subregion  creation  phase.  The  indices  of  this  array  are  the  ID's  of 
subregions.  It  keeps  necessary  information  about  the  subregions  such  as  the  equation 
of  the  plane  fitted  to  the  subregion,  the  boundary*  of  the  subregion,  and  other  subsidiary 
information. 

The  general  structure  of  my  first  two  modeling  algorithms  is  shown  in  Figure  1  on 
page  12. 
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Figure  1.      Block  Diagram  of  the  First  Two  Program  Approaches 

B.     JOINT  TOP-DOWN  AND  BOTTOM-UP  PLANAR-PATCH  TERRAIN 
MODELING 

My  first  algorithm  consists  of  two  parts.  The  first  part  is  the  top-down  subdivision 
of  the  original  elevation  data  matrix  into  subregions  and  the  second  part  is  the 
bottom-up  merging  of  adjoining  subregions  when  adjoining  subregions  have  similar 
plane  equations. 
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1.  TOP-DOWN  PHASE 

A  quadtree  region  subdivision  method  is  used  during  the  top-down  phase  of  the 
firtst  algorithm.  The  method  generates  a  quadtree  by  recursively  dividing  a  two- 
dimensional  region  into  quadrants.  Each  node  in  the  quadtree  has  four  data  elements, 
one  for  each  of  the  quadrants  in  the  region.  [Ref  4:  p.  215]  The  original  region  is  given 
some  standard  tests  such  as  comparison  to  a  threshold  of  the  maximum  magnitude  or 
the  minimum  magnitude  of  the  points  within  it,  depending  on  the  application.  Only  if 
a  region  fails  the  tests  and  is  not  of  minimum  size,  it  is  further  subdivided  into  four  su- 
bregions.  If  a  region  has  both  length  and  width  four  cells  or  more,  it  is  not  of  the  min- 
imum size.  See  Figure  4  on  page  16  for  the  different  kinds  of  minimal  subregions.  The 
subdivision  process  described  above  is  repeated  for  each  subregion.  Figure  2  on  page 
14  shows  graphically  this  algorithm. 

The  code  works  whether  a  region  has  a  dimensions  of  odd  or  even  numbers.  For 
example,  consider  4  by  4  and  5  by  5  matrix  regions  to  subdivide;  Figure  4  on  page  16 
shovv's  how  this  method  works. 

2.  THRESHOLDS  USED 

Three  thresholds  are  used  in  the  top-down  phase:  low_bound,  upper_bound  and 
standard_deviation  (SD). 

The  threshold  upper_bound  is  set  by  the  maximum  climbing  capability  of  the 
vehicle  in  cases  where  the  polyhedral  model  will  be  used  for  route  planning.  If  the  ab- 
solute value  of  the  minimum  slope  of  a  region  is  greater  than  the  upper_bound,  which 
means  every  cell  has  a  slope  greater  than  the  upper_bound.  then  the  region  is  assumed 
being  untravcrsable  by  the  vehicle,  and  the  region  need  not  be  subdivided  further  be- 
cause a  good  polyhedral  fit  just  doesn't  matter  for  that  region.  The  threshold  low_bound 
means  that  when  the  absolute  value  of  the  maximum  slope  of  a  region  is  less  than  the 
low_bound,  which  means  ever>'  cell  has  a  slope  less  than  the  lo>v_bound,  then  the  slopes 
in  the  region  are  unimportant.  Such  regions  need  not  be  subdivided  either. 

For  the  region  which  has  maximum  slope  greater  than  the  lo\v_bound  and  mini- 
mum slope  less  than  the  upper_bound,  the  region  is  subdivided  in  four  parts  according 
to  the  quadtree  algorithm  to  better  isolate  the  unlevel  slopes.  Then  we  fit  a  plane  to  the 
region  and  calculate  the  standard  deviation  of  the  plane  from  the  real  terrain  points. 
The  standard  deviation  (SD)  is: 
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QUADTREE  REGION  SUBDIVISION  METHOD 


AN  EXAMPLE  OF  SUBREGIONS  CREATED  BY  THE  QUADTREE  METHOD 


Figure  2.      Quadtree  Region  Subdivision  Method 
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PROCEDURE  MAKE_SUBREGIOHS(YL,YU,XL,XU  =  INTEGER; 

prev_region    region_ptr); 

(X  YL  :  LOWER  BOUND  IN  Y  AXIS,  YU  ••  UPPER  BOUND  IN  Y  AXIS, 
XL  :  LONER  BOUND  IN  X  AXIS,  XU  =  UPPER  BOUND  IN  X  AXIS  X) 

VAR 

P  :  REGION_PTR; 

BEGIN 

IF  ((YU-YL)  >  2)  AND  ((XU-XL)  >  2)  THEN  (X  NOT  MINIMUM  SIZE  X) 

if  ( p3 . max_mag_of_region  >  low_bound)  then     (x  test  x) 
begin 

nake_subregions( ( (yu-yl)  div  2)+yl+l  ,  yu, 

((xu-xl)  div  2)+xl+l  ,  XU  ,  p);     (x  1st  quadrant  X) 
tnake_subregions(yl  ,  ((yu-yl)  div  2)+yl, 

((xu-xl)  div  2)+xl+l  ,  XU  ,  p);     (x  2nd  quadrant  X) 
make_subregions(yl  ,  ((yu-yl)  div  2)+yl, 

xl  ,  ((xu-xl)  div  2)+xl,  p);       (x  3rd  quadrant  X) 
make_subregions( ( (yu-yl )  div  2)+yl+l  ,  yu, 

xl  ,  ((xu-xl)  div  2)+xl  ,  p);      (x  <4th  quadrant  X) 

p  :=  prev_region  ;  (X  goto  parent  region  x) 

end 

else  stop_subdivide   (X  region  passed  standard  test  X) 
else  stop_subdivide      (X  region  has  miniinum  size      X) 

END 


Figure  3.      Pascal  Code  for  tlie  Quadtree  Region  Subdivision  Method 

where  y(.vj-)  is  the  elevation  of  point  (x,y) 

ax  +  by  ^  c  is  the  plane  equation  of  the  region 
n  is  the  total  number  of  points  in  that  region. 

The  region  is  subdivided  only  if  the  calculated  SD  is  bigger  than  the  threshold 
SD  and  it  does  not  have  the  minimum  size. 
3.     BOTTOM-UP  PHASE 

We  now  merge  adjoining  regions  having  similar  plane  expressions.  The  reason 
is  that  in  the  process  of  applying  the  quadtree  subdivision,  we  often  subdivide  a  region 
which  should  not  have  been  so  completely  subdivided.  For  example,  say  a  region  has 
very  steep  slopes  in  the  first  quadrant  and  the  rest  of  the  quadrants  are  flat.  Since  the 
original  region  will  fail  in  the  maximum  slope  test  due  to  the  slopes  in  the  first  quadrant, 
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DIFFERENT  KINDS  OF  MINIMUM  SIZE 
SUBREGIONS 


4 


BOUNDRRIES 


Figure  4.  Minimum-Size  Subregions  Created  by  the  Quadtree  Method:  This  pic- 
ture shows  quadtree  subdivision  method  apphed  to  the  regions  which 
have  dimensions  of  odd  or  even  numbers,  and  also  shows  boundaries 
and  ditTcrent  kinds  of  nunimum  size  subregions. 

the  top-down  phase  will  subdivide  it  in  four.     However,  the  first,  second  and  third 
quadrants  are  all  flat  and  should  not  be  separate. 

The  adjoining  regions  of  a  region  are  found  from  the  point_array  by  looking  up 
the  subregion  ID  field  of  adjoining  point.  The  bottom-up  phase  is  done  in  two  sub- 
phases:  row^combine  and  col_combine.  The  row_combine  finds  and  merges  adjoining  re- 
gions which  are  located  in  a  same  row,  and  the  coI_comnine  finds  and  merges  the  regions 
in  a  same  column. 
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Figure  5.       Subregions  after  the  Merging  Process 


These  subregions  are  created 
by  performing  the  merging  process  over  the  initial  subregions  created  by 
the  quadtree  algorithm. 


To  merge  adjoining  subregions  which  have  similar  plane  equations,  a  value 
equation_difference  is  calculated  from  the  two  plane  equations  and  compared  to  the 
threshold  abc_difference.    The  equation_difference  is: 


equal  ion  _differefK€  =  \j  [a^  —  02)'  "*"  (^:  ~  ^2)    +  (^^i  ""  ^2)" 

where  a^,  b^  and  Cj  are  the  coefTicients  of  plane^  and 
a,.  b2  and  c,  are  the  coefTicients  o{^ planej. 

If  this  is  less  than  the  threshold  then  the  two  adjoining  regions  pass  the  test. 
Then  we  fit  a  plane  over  the  two  subregions  and  calculate  the  standard  deviation  of  the 
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plane,  and  only  if  it  is  less  than  the  threshold  standard_deviation  we  do  merge  the  two 
subregions. 

C.     STRICT  BOTTOM-UP  PLANAR-PATCH  TERRAIN  MODELING 

1.  POINT-CLASSIFICATION  PHASE 

In  strict  bottom-up  modeling,  we  first  classify  every  point  in  a  region  by  at- 
taching two  tags  to  it.  The  first  tag  classifies  slope  of  the  point:  level,  safe-slope  or 
unsafe-slope.  A  threshold  lo>v_bound  defines  the  low  limit  of  slope  and  another  threshold 
upper_bound  defines  the  upper  limit.  The  second  tag  classifies  the  curvature  of  the  point, 
by  the  terrain  cell  classification  criteria  described  in  chapter  2.  The  designators  are  hill, 
depression,  plane  and  special.  The  tag  special  refers  to  all  the  rest  of  the  cell  classifica- 
tions except  the  above  three.  A  threshold  called  cun^ature  is  compared  to  the  eigenvalue 
of  the  Hessian  matrix;  if  its  absolute  value  is  less  than  the  threshold,  then  the  eigenvalue 
is  regarded  as  zero.[Ref  3:  p.  69]  We  group  the  contuguous  points  with  the  same  first 
and  second  tag  to  form  subregions. 

Next  the  thresholds  magnitude_ratio  and  abc_difference  are  used  to  combine  the 
subregions.  The  magnitude_ratio,  is  used  in  the  region-growing  phase  and  the 
abc_difference  is  used  as  in  joint  top-down  and  mottom-up  modeling. 

2.  REGION-GROWING  PHASE 

The  region-growing  phase  consists  of  two  subphases:  one-dimensional  and 
two-dimensional. 

a.     ONE-DIMENSIONAL  REGION-GROWING 

The  one-dimensional  phase  starts  with  the  first  point  (p^^)  in  the  first  row 
of  the  point_array.  The  point  becomes  the  head  of  a  linked  list  and  the  subregion_ID 
field  of  the  point  is  given  the  value  one.  The  subregion  ID  starts  from  one  and  is  in- 
cremented every  time  a  new  subregion  is  created.  Whenever  a  new  subregion  is  created, 
the  necessary'  information  about  the  subregion,  such  as  the  coordinates  of  the  starting 
point  and  boolean  value  that  says  whether  the  subregion  is  active  or  not,  is  stored  in  the 
subregion_array  using  the  subregion  ID  as  an  index.    Then  it  examines  the  next  point  ( 

/7p)  in  the  same  row.  If  the  ratio  computed  using  magnitudes  of  gradient  ( 
aDs{magnnude^^^  —  magnitude ^xx) 

'■ — ; — ^^ )  is  less  than  the  threshold  magnitude  ratio,  then  they 

niagniiude^-^^ 

are  linked  together  to  form  a  Hnear  subregion  (/^n)-  That  is,  the  next__point  field  of  the 
first  point  in  the  point_array  is  given  the  x  and  y  coordinate  of  the  second  point  and  the 
subregion_ID  field  of  the  second  point  is  given  the  same  value  as  the  subregion  ID  of  the 
first  point.    The  average  magnitude  of  the  gradient  of  the  points  within  the  subregion  is 
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Figure  6.      One-Dimensional  Re"ion-Gro>>inj! 

calculated  and  is  stored  in  the  avg_niag  field  of  the  subregion  record  in  the 
subregion_array.  Then  repeatedly  for  each  next  point  {f\j  in  the  same  row.  the  ratio 
between  the  average  magnitude  of  the  gradient  of  the  subregion  and  the  magnitude  of 
the  gradient  of  the  new  point  is  compared  to  magnitude_ratio.  If  on  this  basis  new  point 
(/'iJ  should  not  join  the  subregion,  a  new  region  starts  from  that  new  point  and  the  new 
point  becomes  a  head  of  another  linked  list.  This  process  continues  for  all  rows.  See 
Figure  6  for  an  example  of  this  algorithm. 

b.     TWO-DIMENSIONAL  REGION-GROWING 

The  two-dimensional  subphase  starts  with  the  first  linear  subregion  (/?,i) 
which  was  created  by  the  one-dimensional  subphase.  Actually,  it  starts  from  the  first 
point  in  the  first  column  of  the  point_array.  which  belongs  to  the  subregion  one.  It  ex- 
amines the  point  below  in  the  same  column,  which  belongs  to  another  subregion,  to 
decide  whether  the  average  magnitude  of  the  gradient  of  the  subregions  are  similar.  The 
information,  average  magnitudes  of  the  gradient  of  the  two  subregions,  is  retrieved  fi^om 
the  the  subregion_array  using  the  subregion  ID  of  each  point.  If  the  ratio  is  less  than 
magnitude_ratio  then  they  are  linked  together.    The  tail  o^  the  first  subregion  is  linked 
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Figure  7.      TA>o-Dimensional  Region-Groning:     This  algorithm  is  apphed  after  the 
one-dimensional  rcgion-growine  algorithm. 

to  the  head  of  the  second  subregion.  the  tail  of  the  second  subregion  is  left  nil,  all  the 
subregion  IDs  of  the  points  within  the  second  subregion  are  changed  to  the  first  subre- 
gion ID,  and  the  second  subregion  in  the  subregion_array  is  marked  inactive.  When 
those  two  linear  subrcgions  are  merged,  the  resulting  subregion  is  not  one-dimensional 
any  more,  except  if  both  subregions  were  single  points.  Then  repeatedly  for  each  next 
point  in  the  same  column,  the  ratio  between  the  average  magnitudes  of  the  gradient  of 
the  subregions  is  compared  to  the  niagnitude_ratio.  If  a  new  subregion  does  not  join  the 
old  subregion.  the  control  is  moved  to  the  new  subregion  and  the  process  starts  from 
that  new  subregion.  This  process  continues  for  all  columns.  This  is  illustrated  in 
Figure  7.     After  the  two-dimensional  phase,  we  fit  a  plane  to  each  subregion  using 
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least-squares.  Then  the  threshold  abc_difference  is  used  to  merge  adjoining  regions  that 
have  similar  plane  expressions.  The  merging  is  done  in  two  steps.  The  first  step  is  to 
scan  each  row  of  the  point_array  looking  for  adjacent  points  in  the  same  row  which  have 
different  subregion  IDs.  For  ever\'  pair  of  adjacent  subregions,  the  plane  expressions 
are  consulted  from  the  subregion_array  by  using  the  ID  numbers  of  the  subregions  as 
the  indices.  The  equation_difference  of  the  two  plane  equations  are  calculated  and 
compared  to  abc_difference.  If  it  less,  the  two  subregions  are  merged,  point  of  the  sec- 
ond subregion  to  the  tail  point  of  the  first  subregion.  The  second  step  is  identical  except 
that  scanning  for  merging  is  done  by  column. 

D.     AN  ALGORITHM  IMPROVING  THE  CONTINUITY  OF  ADJOINING 
REGIONS 

One  common  problem  with  planar-patch  terrain  modeling  is  that  planar  patches  for 
two  adjacent  subregions  seldom  yield  a  surface  that  is  continuous  along  the  line  seg- 
ments bisecting  the  points  used  to  create  the  regions.  Since  the  initial  requirements  for 
this  research  rule  out  any  other  terrain  modeling  than  planar-patch,  we  had  to  come  up 
with  some  idea  to  ensure  continuity  of  the  patches.  One  simple  solution  is  to  model 
terrain  entirely  with  small  triangles,  each  of  which  exactly  fits  three  adjacent  terrain 
points,  as  described  in  Chapter  1.   The  number  of  triangles  necessary  is  very  large. 

But  a  minor  modification  to  the  two  modeling  approaches  explored  in  this  research 
can  give  smoothness.  One  can  create  supplementary  triangular  subregions  to  fill  the  gap 
between  adjoining  subregions.  The  gap  between  adjoining  subregions  have  two  lines,  one 
for  each  side  of  the  gap.  The  lines  are  the  edges  of  the  adjacent  planes.  We  can  define 
two  triangles  on  the  gap  by  grouping  the  first  end  point  of  the  first  edge  to  the  second 
edge  and  second  point  of  the  second  edge  to  the  first  edge.   See  Figure  8  on  page  22. 
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TWO  TRIANGLES  CONNECTING 
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SUBREGION  1       SUBREGION  2 


Figure  8.      Fillinjz  the  Boundan  Avith  Two  Triangles 


This  is  better  than  the  many-triangle  model,  but  is  still  not  optimal,  because  the  number 
of  subregions  v\ill  be  increased  considerably  by  the  supplementary  triangular  subregions. 

We  have  made  a  difTcrent  attempt  to  solve  this  problem  using  smoothing  of  the  el- 
evation array  and  an  expansion  to  three  dimensions  of  Rolle's  theorem.  The  basic  idea 
is  that  if  we  fit  planar-patches  to  an  ever\"\vhere-convcx  or  ever>'\vhere-concave  region, 
the  fitted  patches  are  much  more  likely  to  be  continuous.  This  means  we  should  try  to 
ensure  that  all  the  points  within  large  regions  are  classified  as  either  planar  points,  hill 
points  or  depression  points.  See  Figure  9  on  page  23.  Smoothing  t-he  elevation  data 
will  help:    other  kinds  of  points  will  tend  to  be  changed  into  one  of  the  three. 

Unfortunately,  because  of  time,  we  were  not  able  to  solve  this  rather  complex 
problem  completely.    It  is  recommended  for  future  research. 
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INTERSECTION 


THE  INTERSECTION  TENDS  TO  LIE  BETWEEN  THE  TWO  REGIONS 
IF  WE  FIT  PLANES  OVER  THE  EVERY-WHERE  CONVEX  OR 
EVERY-WHERE  CONCAVE  TERRAIN 


Figure  9.      Ho^^  Surface  Convexity  Helps 
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IV.     DISCUSSION 

A.  ACCURACY  ANALYSIS 

The  joint  top-down  and  bottom-up  and  strict  bottom-up  approaches  produce  fairly 
good  terrain  models  in  terms  of  accuracy  and  simplicity.  Here  accuracy  means  whether 
the  model  represented  the  real  terrain  without  losing  anything  significant,  and  simplicity 
means  the  number  of  subregions. 

The  overall  accuracy  of  a  model  can  be  verified  by  comparing  the  pictures  of  the  real 
terrain  and  model.  But  there  is  another  thing  that  describes  the  accuracy  of  a  model, 
standard  deviation.  We  can  evaluate  the  local  accuracy  in  terms  of  the  standard  devi- 
ation of  the  elevation  dilTerence  for  a  point  in  the  real  terrain  and  the  corresponding 
point  in  the  model.  In  the  top-down  approach,  each  subregion  must  have  a  fit  standard 
deviation  less  than  the  the  threshold  standard_deviation.  Let's  assume  that  the  elevation 
difierence  between  a  point  on  the  real  terrain  and  a  point  on  the  model  has  a  normal 
probability  distribution.  If  we  know  the  standard  deviation  of  the  elevation  differences 
of  a  subregion: 

1.  About  68  percent  of  the  points  will  have  elevation  difference  less  than  plus  or  mi- 
nus one  standard  deviation. 

2.  About  98  percent  of  the  points  will  have  elevation  difference  less  than  plus  or  mi- 
nus two  standard  deviations. 

3.  Practically  all  points  (99.73  percent) 

will  have  elevation  difference  less  than  plus  or  minus  three  standard  deviations. 

So  for  example,  if  we  limit  the  model  to  have  the  worst  case  local  elevation  differ- 
ence, say  10  feet,  with  98  percent  certainty,  we  can  set  the  standard  deviation  threshold 
to  5  feet.  Consequently  everv'  subregion  will  have  a  standard  deviation  less  than  5  feet, 
and  98  percent  of  the  points  will  have  an  elevation  difference  less  than  10  feet. 

The  threshold  standard  deviation  was  not  used  in  the  strict  bottom-up  approach  so 
the  above  argument  does  not  hold. 

B.  LIMITATIONS  OF  THE  MODEL 

In  the  joint  top-down  and  bottom-up  approach,  the  minimum  number  of  points  in 
a  subregion  is  four  and  the  shape  of  initial  subregion  must  be  a  rectangle.  This  re- 
striction is  due  to  the  quadtree.  These  limitations  prohibit  the  approach  from  picking 
up  linear  (one-cell- wide)  terrain  features  and  creating  linear  subregions.   This  limitation. 
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however,  does  not  exist  in  the  strict  bottom-up  approach  owing  to  the  hnked  list  of 
points  used  to  represent  a  region.  When  a  hnear  region  is  a  straight  horizontal  or  ver- 
tical line,  one  cannot  calculate  a  plane  equation. 

C.     EXPERIMENTAL  RESULTS 

We  now  present  the  output  of  our  terrain  modeling  program  in  graphical  form.  It  is 
apphed  to  sample  elevation  data  acquired  from  the  Fort  Hunter  Liggett  terrain  data 
base.  Shown  are  the  pictures  that  are  produced  by  our  modeling  programs  using  differ- 
ent thresholds. 

As  was  mentioned  in  chapter  3,  in  joint  top-down  and  bottom-up  modeling,  the 
thresholds  were  Io>v_bound  and  upper__bound  for  gradient  magnitude,  standard_deviation, 
and  abc_difference.  In  strict  bottom-up  modeling,  the  thresholds  were  lo>v_bound  and 
upper_bound  of  gradient  magnitude,  curvature,  abc_difference,  and  magnitude_ratio. 

Three  sets  of  pictures  are  given  sequentially  in  the  following  section.  The  first  set 
is  produced  by  the  joint  top-down  and  bottom-up  approach  (referred  as  the  first  ap- 
proach in  the  following  pictures),  the  next  by  the  strict  bottom-up  approach  (referred 
as  the  second  approach),  and  the  last  by  the  bottom-up  approach  applied  to  the 
smoothed  data  (referred  as  the  third  approach). 

In  the  array  pictures,  each  subregion  has  a  unique  code.  If  a  number  does  not  have 
a  prefix  then  the  number  is  actually  the  subregion  ID.  If  it  does  have  a  prefix  then  the 
number  has  to  be  converted.  See  the  conversion  example  below. 

23  =  =  >  subregion  ID   23. 

.23  =  =  >  subregion  ID  123. 

:23  =  =  >  subregion  ID  223. 

-23  =  =  >  subregion  ID  323. 

0  =  =  >  either  an  outermost  boundary 

or  inside  area  of  a  region. 
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1.     EXPERIMENTAL  RESULTS  1  OF  THE  FIRST  APPROACH 
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Figure  10.      Sample  Real  Terrain:     This  is  the  real  terrain  that  were  used  in  these 
experiments. 
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Figure   II. 


Subregions     Created     by     the     Quadtree     Method     in     the     First 
Approach:     Thresholds:  lo\v_bound  =  0.1;       upper_bound  =  0.6; 

siandard_deviation=  3.  The  above  picture  shows  the  initial  subregions 
created  by  quadtree  algorithm. 
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Figure  12.  Subregions  Created  after  the  Merging  in  the  First  Approach:  Thresh- 
olds: abc_difrerence=  10;  standard_deviation  =  3;  Some  of  the  initial 
subregions  created  by  quadtree  algorithm  are  merged  by  merging 
process,  thus  yielding  less  subregions. 
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Figure   13.      Terrain  Model  Created  by  the  First  Approach:      This  is  the  final  prod- 
uct of  the  first  approach. 
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Figure  14.      Terrain  Model  Created  by  the  First  Approach  Showing  Boundaries 
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Subregions     Created     by     the     Quadtree     Method     in     the     First 
Approach:     Thresholds:  low_bound  =  O.I;       upper_hound  =  0.6; 

standard_dcviation=  7.  The  above  picture  shows  imiial  subregions  cre- 
ated by  the  quadtree  algorithm  using  a  different  threshold  (the 
standard_deviation  is  chansed). 
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Figure  17.      Terrain  Model  Created  by  the  First  Approach:      Final  product  of  the 
first  approach  using  a  diilerent  set  of  thresholds. 
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Figure  18.  Terrain  Model  Created  by  the  First  Approach  Showing 
Boundaries:  The  planar  patches  are  difTerent  from  the  planar  patches 
produced  from  the  previous  experiment. 
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Figure  19.      Subregions  Created  by  the  Region  Growing  Algorithm  in  the  Second 
Approach:     Thresholds:  low_bound  =  0.1;       uppcr_bound  =  0.6; 

magmtude_ratio  =  0.1.    This  picture  is  produced  after  performing  the 
one-dimensional  region  growing  process  (row  region  growing). 
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Approach:  Threshold:  abc_di(Terence=  10;  the  two-dimensional  re- 
gion growing  algorithm  is  applied  over  the  subregions  created  by  the 
one-dimensional  region  growing  algorithm  and  then  merging  is  per- 
formed. 
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Figure  21.      Terrain  Model  Created  by  the  Second  Approach:     Final  product  of  the 
second  approach. 
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Figure  22.      Terrain  Model  Created  by  the  Second  Approach  Showing  Boundaries 
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■4.     EXPERIMENTAL  RESULTS  2  Of  THE  SECOND  APPROACH 
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Figure  25.  Terrain  Model  Created  by  the  Second  Approach:  Final  product  of  the 
second  approach.  Notice  that  this  picture  is  diflerent  from  the  picture 
produced  by  the  previous  experiment. 
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Figure  26.      Terrain  Model  Created  by  the  Second  Apprpach  Showing  Boundaries 
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5.     EXPERIMENTAL  RESULTS  OF  THE  THIRD  APPROACH 
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Figure  27.      Smoothed  Sample  Terrain:     The  original  sample  terrain  is  smoothed. 
Bottom-up  modeling  vvill  be  performed  on  this  smoothed  terrain. 
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Figure  28.  Subregions      Created      by      the      Third      Approach:     Thresholds: 

lo\v_bound  =  0.1:  upper_bound  =  0.6;  magniiude_raiio  =  0.1. 
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Figure  29.  Subregions      Created      after      the      Merging      in      the      Third 

Approach:     Threshold:   abc_difrercnce  =  10. 
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Figure  30.  Terrain  Model  Created  by  the  Third  Approach:  This  picture  looks 
ahnost  the  same  as  the  the  picture  produced  by  the  second  approach; 
this  approach  produced  less  subregions  than  the  second  approach. 
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Figure  31.      Terrain  Model  Created  by  the  Third  Approach  Showing  Boundaries 
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V.     CONCLUSIONS 

In  this  research,  we  explored  three  types  of  planar-patch  terrain  modeling.  The  purpose 
of  such  models  is  a  simple  and  yet  accurate  view  of  terrain.  Since  our  models  were  de- 
veloped for  route  planning,  we  required  a  minimal  number  of  subregions,  since  this  di- 
rectly affects  the  cost  of  the  route  planning. 

The  first  algorithm  (joint  top-down  and  bottom- up  algorithm)  takes  the  longest  ex- 
ecution time  among  the  three.  The  second  (strict  bottom-up)  and  the  third  (bottom-up 
algorithm  applied  to  the  smoothed  elevation  data)  algorithms  take  almost  same  time, 
though  the  third  algorithm  requires  the  input  elevation  data  to  be  smoothed  first,  re- 
quiring some  extra  preprocessing  time.  As  far  as  result  accuracy  (simplicity)  is  con- 
cerned, the  first  algorithm  is  better  than  the  others.  The  second  and  third  do  not  use  a 
standard  deviation  threshold,  resulting  in  possible  inaccuracy,  and  tend  to  produce  more 
subregions.  Since  all  of  the  three  algorithms  use  almost  the  same  data  structure,  a  two- 
dimensional  point_array  and  a  one-dimensional  subregion_array,  they  take  almost  the 
same  space.  As  for  the  application  environments  of  these  algorithms,  the  first  algorithm 
could  be  used  in  situations  where  the  accuracy  is  critical,  while  the  second  and  third 
algoritms  could  be  preferred  in  situations  where  a  short  execution  time  is  required  and 
linear  terrain  features  should  be  identified. 

Since  the  algorithms  used  in  this  research  operate  on  large  arrays  and  involve  hea\'\' 
computations,  a  multiple-processor  computer  architecture  would  yield  superior  process- 
ing capability  since  the  same  operations  could  be  performed  for  each  part  of  the  input 
matrix  simumtaneously.  As  I  conclude  this  research,  I  regret  that  I  could  not  solve  the 
continuity  issue  completely. 
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APPENDIX  A.     PASCAL  CODE  OF  JOINT  TOP-DOWN  AND 
BOTTOM-UP  ALGORITHM 

(*$s350000''^) 

Author   =  Yee,  Seung  Hee. 
Date     =  3  June  1988. 

Input    =  1.  Elevation  data  acquired  from  Fort.  Hunter  Liggett 
data  base  (40  by  40  square). 
2.  file  of  magnitude  of  gradient  of  the  above  elevation 
data  created  by  using  the  program  in  Appendix  3. 
Output   =  1.  Elevation  data  of  the  model  created  by  this  program. 
2.  Elevation  data  of  the  boundary  points,  other  points 
are  given  the  value  one  so  that  when  we  display 
this  model  in  three-dimensional  pictures,  we  can 
only  see  the  boundaries  of  the  patches. 
Computer  =  Waterloo  Pascal  on  IBM  370/3033AP, 

Naval  Postgraduate  School. 
Description  = 

This  program  is  the  Pascal  implementation  of  the 
Joint  top-down  and  bottom-up  terrain  modeling. 


program  top_down( input, output); 
const 

mtinft  =  3.  2808; 

low_bound  =  0.  1;    ('"^  mag  >  10  %   slope   then  divide  *) 
{">''   mag  <=  10  %   then  disregard  *) 

upper_bound  =  0. 6;    (*  mag  <  60  %   then  do  not  subdivide  *) 

max_sd   =  3; 

abc_dif ference  =10  ; 

max_group_nura=  300; 

row  =  40; 

col  =  40; 

increment  =  3; 
type 

elrange  =  1. . 10000; 
pointrec  =  record 


el  : 

e 

1  range  ; 

mag: 

real; 

gP  ' 

0 

. max_group_num; 

end; 

magrec  = 

record 

mag 

; 

real; 

r,c 

: 

1.  .  row; 

end; 

point_array=array  (.  1.  .  row,  1.  .  col.  )  of  pointrec; 
grouppointer  =  -'grouprec  ; 
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grouprec  =  record 

r 1 , rh , 

cl  ,cr  :  1.  .  row; 

magmax,magmin  :  magrec  ; 

up    :  grouppointer; 
end; 
constarray  =  array  (.  1.  .  3, 1.  .  4.  )  of  real; 

nodepointer  =  --node; 
node  =  record 

r ,c  : integer; 

next  :  nodepointer; 
end; 

coefrec  =  record 

crljCrh, 

ccljccr  :  1.  . row; 
a,b,c  :  real  ; 
minmag,maxmag  :  real; 
sd   :  real  ; 
ee,ex,ey  :  real; 
side  :  0. . max_group_num; 
constarr  :  constarray; 
active  :  boolean; 
list   :  nodepointer  ; 
end; 

subregion_array  =  array  (.  1. . max_group_num.  )  of  coefrec  ; 


var 

points  :  point_array  ; 

groups  :  grouprec  ; 

subregion  :  subregion_array  ; 

ok     :  constarray; 

data,datara,doutl ,dout2  :  text; 

counter ,r,c, i  :  integer; 

avg,ee,ex,ey   :  real; 

top,p  :  grouppointer; 

temparr:  array  ( . 1. . max_group_num. )  of  boolean; 

function  variance  (a,b,c,ee,ex,ey: real;  ok: constarray): real; 

begin 

variance  :=  (ee-(2*a*ex)-(2*b*ey)-(2*c*ok(.  3,4.  )) 
+( a^^a^>-ok( .  1 , 1.  ) )+( 2''-'a'^b''^ok(  .2,1.)) 
+( 2'«b^^c-'^ok( .  2 , 3.  )  )+(  2*a--c^-ok(  .1,3.)) 
+(b^--b''-ok(.2,2.  ))+(c*c*ok(.3,3.  ))) 
/(ok(.3,3.)-l); 

end; 

procedure  initialize; 
var 

r,c  : integer; 
begin 

counter  :  =  0  ; 

for  c  :  =  1  to  col  do 
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for  r  : =  1  to  row  do 

points(.  r.  )(.  c.  )•  gP  :=  0  ; 
for  r  : =  1  to  max_group_num  do 
begin 

subregion( .  r.  )•  a  :=  -999; 

subregion( .  r.  )•  b  :=  -999; 

subregion( .  r.  )•  c  :=  -999; 

subregion( .  r.  ).  active  :=  false; 

subregion( .  r.  )•  side  :=  0; 

subregion(.  r.  )•  sd  :=  -999; 


end; 


end; 


(-v  procedure  initialize  *) 


procedure  getconst(var  points: point_array;  rl ,rh,cl ,cr: integer; 

var  arr:  constarray;  var  ee,ex,ey  : real); 
var 

i, j:  integer; 


begin 
for 


i  : =  1  to  3  do 
begin 

for  j  :=  1  to  4  do 
arr(.  i,j.  ):=0; 
end; 
ee  :=  0; 
ex  :=  0; 
ey  :  =  0; 

for  i  : =  rl  to  rh  do 
begin 

for  j  : =  cl  to  cr  do 
begin 


arr(.  1,1.) 
arr(.l,2.  ) 
arr(.  1,3.  ) 
arr(. 1,4. ) 


=arr(.  1,1.  )+(i'U); 

=arr(.l,2.  )+(j*i); 

=arr(.  1,3.  )+  i; 

=arr(.  1,4.  )+(points(.  i,  j.  ).  el'^^i); 


arr(.2,l.  ) 
arr(.2,2.  ) 
arr(.2,3.  ) 
arr(.2,4.  ) 


=arr(.2,l.  )+(i*j); 

=arr(.2,2.  )+(j''^j); 

=arr(.2,3.  )+  j; 

=arr(.  2,4.  )+(points(.  i,  j.  ).  61'"^^); 


arr(. 3,1.  ) 
arr(. 3,2.  ) 
arr(. 3,3.  ) 
arr(.3,4.  ) 


=arr(.  3,1.  )+i; 
=arr(.3,2.  )+j; 
=arr(.  3,3.  )  +  l; 
=arr(.  3,4.  )+points(, i,j. ). el; 


ee  :=  ee  +  sqr(points( .  i,  j.  ).  el); 
ex  :=  ex  +  ( points (.  i,j.  ).  el  *  i  ); 
ey  :=  ey  +  (points( .  i,  j.  ).  el  *  j  ); 

end  (*  j  *) 
end;       (^'  i  ^0 
end;  (,*   procedure  get  const  *) 
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procedure  gauss(ok: constarray;  var  aa,bb,cc: real); 
var 


i,j  :  integer; 

temp  :  array  (.1..4.  )  of  real; 
t  :  real; 
begin 


(*  error  checking  purpose  *) 


i  :=  1; 

while  ok(.i,l.  )  =  0.  0  do  i  :=  i  +  1  ;  (*  row  rotation  *) 

for  j  :=  1  to  4  do 

begin 


temp(. j.  ) 
ok(.l,j.  ) 
ok(. i,j.  ) 

end; 

t  :=  ok(.l,4.  ); 


=  ok(.l,j.) 
=  ok(. i,j. ) 
=  terapC. j.  ) 


for  i  :  =  2  to  3  do 

for  j  :  =  4  downto  1  do 

ok(.i,j.)  :=  ok(.i,j.)'^ok(.l,l.  )  +  ok(.l,j.)*(-ok(.i,l.  )); 


i  :=  2; 

while  ok(.i,2.)  =  0.  0  do   i 

for  j  :=  1  to  4  do 

begin 


=  i  +  1  ; 


temp(. j.  ) 
ok(.2,j.) 
ok(.  ij.  ) 


=  ok(.2,j.) 

=  ok(.  i,j.  ) 
=  temp(. j.  ) 


end; 


for  j  : =  4  downto  1  do 

ok(.3,j.)  :=  ok(.3,j.)'^ok(.2,2.  )  +  ok( .  2  ,  j.  )*( -ok( .  3,2.  ) ); 

if  ok(. 3,3. )  <>  0. 0  then 
begin 


ok(.  3,4. ) 
ok(.2,4.  ) 
ok(.l,4.  ) 


=  okC.  3,4.  )  /  ok(.3,3.  ); 

=  (  ok(.2,4.  )  -  ok(.3,4.  )^^ok(.2,3.  ))  /  ok(.2,2.  ); 

=  (  ok(.l,4.  )  -  ok(.2,4.  )^'ok(.  1,2.  ) 

-  ok(.3,4.  )''^ok(.l,3.  ))  /  ok(.  1,1.)  ; 


("  solution  *) 

("  solution  ") 

('-  solution  ") 
(*  gauss  error  =  t-(   aa^>-ok( .  1, 1.  )+bb*ok(.  1,2.  )+  cc*ok(.  1,3.  ))  *) 
end; 
end;     (*  procedure  gauss  *) 


aa  :  =  ok(.  1,4.) 
bb  :=  ok(.2,4.  ) 
cc  :  =  ok( . 3,4. ) 


procedure  prefix(nurn:  integer); 
begin 

case  ((num  div  100)  mod  10)  of 


write( 
write( 
write( 
write( 
write( 
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5 

write( 

%') 

6 

write( 

,V) 

7 

write( 

&•) 

8 

writeC 

?') 

9 

write( 

=  ') 

end; 
end;      (*  procedure  prefix  ''0 

procedure  printsubl(points: point_array); 
var 

r,c,i,j,k  :  integer; 
begin 
page; 

for  k  :  =  1  to  counter  do 
with  subregion( . k. )  do 

for  i  :=  crl+1  to  crh-1  do 
for  j  :=  ccl+1  to  ccr-1  do 
points(.  i,  j.  )•  gp  :=  0; 
for  r  : =  1  to  row  do 
for  c  :  =  1  to  col   do 
begin 

if  c<>  col   then  begin 

prefix(points(.  r,c.  ).  gp); 
write(  (points(.  r,c.  )•  gp  naod   100):  2); 
end 
else  begin 

prefix( points ( . r,c.  ).  gp); 
writeln( (points( . r ,c.  )•  gp  rood  100): 2); 
writeln; 
end; 
end; 
end;      (*   procedure  printsubl  *) 


procedure  printsub2(points: point_array); 
var 

k, r , c,rl , r2 ,cl ,c2  :  integer; 

current  :  nodepointer; 
begin 

for  r:  =  2  to  row-1  do 

for  c:  =  2  to  col-1  do 

point s (.  r,c.  )•  gp  :=  0  ; 

for  k:  =  1  to  counter  do 

if  subregion(. k. )• active  then 

begin 

current  :=  subregion( . k. )• list  ; 
while  current  <>  nil  do  begin 
rl:  =current-'.  r; 
cl:  =current-'.  c; 

if  current--,  next  <>  nil  then  begin 
r2:  =current->.  next--,  r; 
c2:  =current-'.  next--,  c; 
end 

else  begin 
r2  :=  rl; 
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c2  :=  cl; 

end; 

if  rl  >  r2  then  begin  r:  =  r2;  r2:=  rl;  rl:=  r;  end; 


if  cl  >  c2  then  begin  c:  =  c2;  c2:  =  cl;  cl:=  c;  end; 

if  rl=r2  then  for  c  :=  cl  to  c2  do 
points (.  rl,c.  )•  gp  •  =  k; 
if  cl=c2  then  for  r  :=  rl  to  r2  do 
points(.  r,cl.  )•  gp  •' =  k; 
current  :=  current--,  next; 
end; 
end; 

for  r  : =  1  to  row  do 
for  c  : =  1  to  col  do 
begin 

if  c<>  col   then  begin 

pref ix(points(.  r,c.  ).  gp); 
write( (points( . r ,c.  )•  gp  rood  100): 2); 
end 
else  begin 

pref  ix( point s(  .  r,c.  ).  gp); 
writeln( (points( .  r jC.  )•  gp  rood  100): 2); 
writeln; 
end; 
end; 
end; 

(*  procedure  printsub2  *) 


procedure  get_nodes( var  subregion:  subregion_array); 
type 

direction  =  (  up, down, right , left  ); 
var 

i,j,k  :  integer; 

current , head, p  :  nodepointer; 

temp  :  direction  ; 

procedure  get_next(var  p  :  nodepointer; 

var  temp  :  direction); 
var 

i,j  :  integer; 

done  :  boolean  ; 


procedure  do_dir(var  done: boolean  ;  dir:  direction  ); 
begin 

new(p); 


P-.  r 
p--.  c 
done 


=  1  ; 

=  j  ; 

=  true; 


case  dir  of 

up 

down 

right 


temp  :  =  up; 
temp  : =  down; 
temp  : =  right; 
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left   :   temp  :=  left  ; 
end;  (^''  case  '^^ 
end;  {'■'   sub_sub_  procedure  do_dir  *) 


begin  (''<■  sub  procedure  get  next  ''0 
done  :=  false  ; 
i  :=  p-.  r; 
j  :=  P"-  c; 
case  temp  of 

right  :  begin 

while  not  done  do  begin 

if  points(. i-1, j.  ).  gp  =  k 
then  do_dir(done,up); 
if  (not  done  and 

(points(.  i,j  +  l.  ).gp  <>  k)) 
then  do_dir( done, down); 
if(  (  subregion( .  k.  )•  list--,  r  =  i) 

and  (subregion( .  k.  )•  list--,  c  =  j))  then 
begin 

new(p); 
P-.  r  :=  i  ; 
P"-  c  :=  j  ; 
done  : =  true; 
end; 

if  not  done  then  j:=  j+1; 
end;  {*   while  *) 
end; 
left   :  begin 

j  :=  j-1; 

while  not  done  do  begin 

if  points(. i+1, j.  ).  gp  =  k 
then  do_dir( done, down); 
if  (not  done  and 

(points(. i, j-1. )•  gp  <>  k)) 
then  do_dir(done,up); 
if  not  done  then  j:=  j-1; 
end;  ("'^  while  *) 
end; 
up    :  begin 

i  :=  i-1; 

while  not  done  do  begin 

if  points(.  i,  j-1.  ).  gp  =  k 
then  do_dir(done, left); 
if  (not  done  and 

(points(. i-l,j. )-gP  <>  k)) 
then  do_dir( done, right); 
if  not  done  then  i:  =  i-1; 
end;  (^-  while  *) 
end; 
down  :  begin 

i  :=  i+1; 

while  not  done  do  begin 

if  points(.  i,  j+1.  )■  gP  =  k 
then  do_dir( done, right); 


55 


if  (not  done  and 

(points(.  i+l,j.  ).gp  <>  k)) 
then  do_dir( done, left); 
if  not  done  then  i:  =  i+1; 
end;  (■''"  while  *) 
end; 
end;  (*  case  *) 
end;  (*  sub  procedure  get_next  ''0 


begin   (*  procedure  get  nodes-) 
for  k  : =  1  to  counter  do 
if  subregion(. k. )•  active  then 
begin 

new(p); 

p-".  r  :  =subregion(.  k.  )•  crl  ; 

p-".  c  :  =subregion( .  k.  )•  ccl  ; 

while  (  points( .  p-'.  r-ljp-".  c.  )•  gP  -   k)  do 

p--.  r  :=  p--.  r-1;    ("  move  to  the  upper  most  row  of  the  gp  '='0 
while  (  points(.  p--.  rjp--.  c-1.  )•  gp  -   k)  do 

p--.  c  :=  p-".  c-1;    ("  move  to  the  left  most  col  of  the  gp  -) 

subregion(.  k.  )•  list  :=  p  ; 
current  : =  p  ; 
temp  :  =  right  ; 
repeat 

get_next(p , temp) ; 

current--,  next  :=  p  ; 

current  :  =  p  ; 
until  (subregion(.  k.  )-list-'.r  =  p-'.r) 

and  (  subregion( .  k.  )•  list--,  c  =  p-*.  c)  ; 
p--.  next  :  =  nil; 
end; 


(*   for  k  : =  1  to  counter  do 
if  subregion( .  k.  )•  active  then 
begin 

writeCgp  ',k:3,''V.v=>'); 
current  :=  subregion( .  k.  )• list  ; 
while  current  <>  nil  do  begin 

wr it e( current--,  r:  2  ,  '  ,  '  , current--,  c:  2,  '=>'  ); 
current  :=  current--,  next; 
end; 

writeln; 
end*) 
end;  (''''  procedure  get_nodes  *) 


procedure  makesub(rl,rh,cl ,cr  :  integer;  prev: grouppo inter); 
type 

why  =  (nondivisible, magnitude, std, steep); 

var 

dif,r,c  :  integer; 

flag    :  why; 
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sd,a,b,cc  : real; 


(* 


procedure  stopdivide( f lag  :  why  ); 
var 

r,c  :  integer; 
begin 

counter  : =  counter  +  1; 
with  subregion(. counter. )  do 
begin 

crl  :  =  rl; 
crh  :  =  rh; 
ccl  :  =  cl; 
ccr  :  =  cr; 

maxmag  :  =  p-".  magmax.  mag  ; 
minmag  :=  p-".  magrain.  mag  ; 
end; 
for  c  :=  cl  to  cr  do 
for  r  : =  rl  to  rh  do 

points( .  r.  )( •  c.  )•  gp  '•-   counter; 
case  flag  of 

nondivisible  : 

write('"  non  divisible  area!'); 
magnitude 

write('"  stop  divide!  by  mag.  '); 
std  : 

write('"  stop,  it  fits  well.'); 
steep        : 

write('*  stop,  it  is  too  steep.'); 


*) 


end; 

writeln('"  area  rl ,rh,cl ,cr=' ,rl: 3,rh: 3,cl: 3,cr: 3, 
has  gpnum=' , counter: 3  ); 

p  :=  prev  ; 
end;   (-   sub  procedure  stop  divide  *) 


begin    ("  procedure  makesub  *) 
new(p); 

p--.  magmax.  mag  :=  -999; 
p--.  magmin.  mag  :=  999; 


P"- 

rl 

:=  rl; 

p-^- 

rh 

:=  rh; 

P"- 

cl 

:=  cl; 

p"- 

cr 

:=  cr; 

P" 

up 

:=  prev  ; 

for  c  :=  cl  to  cr  do 
for  r  : =  rl  to  rh  do 
begin 

if  points(.  r.  )( .  c.  ).  mag  >  p--.  magmax.  mag  then 
begin 

p-".  magmax.  mag  :  =  points(.  r.  )(.  c.  ).  mag; 
p-^.  magmax.  r  :  =  r; 
p-".  magmax.  c   :  =  c; 
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end; 

if  points(  .  r.  )(  •  c.  )•  fi'Qg  '^  p--.  magmin.  mag  then 

begin 

p--.  magrnin.  mag  :=  points( .  r.  )( •  c.  )•  m^g; 
p--.  magmin.  r   :  =  r; 
p-'.  magmin.  c   :  =  c; 
end; 
end; 
i*       max,min  magnitude  for  region  rl,rh,cl,cr  before  division  *) 

getconstC  points , r 1 , rh , cl , cr , ok , ee , ex , ey ) ; 

gaus  s(ok,a,b,cc); 

sd       :=  sqrt(variance(a,b,cc,ee,ex,ey,ok)); 

(*  sd  for  gp  before  division  *) 

if  ((rh  -  rl)  >  2)  and  ((cr  -  cl)  >  2)  then 
if  (p--.  magmax.  mag  >  low_bound)  and 
(p--.  magmin.  mag  <  upper_bound)  then 
if  sd  >  max_sd  then 
begin 

makesub(rl  ,  ((rh-rl)  div  2)  +  rl, 

cl  ,  ((cr-cl)  div  2)  +  cl,  p); 
makesub((( rh-rl)  div  2)  +  rl  +  1  ,  rh, 
cl  ,  ((cr-cl)  div  2)  +  cl  ,  p); 
makesub(rl  ,  ((rh-rl)  div  2)  +  rl, 

((cr-cl)  div  2)  +  cl  +  1  ,  cr  ,  p); 
makesub((( rh-rl)  div  2)  +  rl  +  1  ,  rh, 

((cr-cl)  div  2)  +  cl  +  1  ,  cr  ,  p); 
p  :=  prev  ; 
end 

else  stopdivide(std) 
else  stopdivide(magnitude) 
else  stopdivide(nondivisible);  (*   smallest  area  *) 
end;  (*  procedure  makesub  *) 

procedure  do_combine(var  pt: point_array; 

var  subregion: subregion_array; 

var  last, new  :  integer); 
var  ii,jj  :  integer; 

tail, current   :  0. . max_group_num; 
begin 

subregion( . last.  ).  active  :=  true; 
subregion( .  new.  ).  active  :=  false; 
current  :=  subregion( . last. ). side  ; 
if  current  =  0  then 
begin 

subregion( .  last.  ).  side  :=  new; 
current  :  =  new; 
end 
else  begin 

while  current  <>  0  do 
begin 

tail  : =  current  ; 

current  :=  subregion(. current.  )•  side  ; 

(*   traverse  *) 
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).  ee 

+ 

).  ee 

).ex 

+ 

).  ex 

)•  ey 

+ 

).ey 

end;  (*   subregion  arr  *) 

subregion(.  tail.  )•  side  :=  new; 

(''f  to  link  new  gp  ''=•) 
current  :=  new;  i*   into  existing  link  *) 

end; 
while  current  <>  0   do  (*  reassigning  gp  #  to  new  region  ''0 
begin 

for  ii  :=  subregion( . current. )• crl  to 
subregion(. current.  ).  crh  do 
for  jj  :=  subregion( . current. ). ccl  to 
subregion( . current. ). ccr  do 
pt(.  iijjj.  ).  gp  :=  last; 
current  :=  subregion(. current. ), side; 
end;   (''f  while  *) 
for  ii  :=  1  to  3  do 
for  j j  : =  1  to  4  do 

subregion( .  last.  ).  cons t arr (.  ii, j j.  )  :  = 

subregion(. last. ). constarr(.  ii,jj. )  + 
subregion(. new. ). constarr( . ii, j j. ); 
subregion( .  last.  ).  ee  :=  coefs(.last. 

subregion( .  new. 
subregion(.  last.  ).  ex  :=  coefs(.last. 

subregion( . new. 
subregion( . last.  ).  ey  :=  coefs(.last. 

subregion( . new. 
with  subregionC . last. )  do 
begin 

gauss (const arr , a,b,c); 

sd       :=  sqrt( variance(a,b,c,ee,ex,ey ,constarr) ); 
end; 
end;     {*   sub  procedure  do_corabine  '''■) 


procedure  combine(var  pt:  point_array; 

var  subregion: subregion_array); 
var 

tabc_difference, i,j , last , new  :  integer; 

joinabcjjoinsd, joined  :  boolean; 

procedure  compareabc(var  join: boolean); 
var 

temp  :  real  ; 
begin 

temp  :=  sqrt(  sqr(subregion( .  new.  ).  a  - 

subregion(. last. ). a)  + 
sqr(subregion( . new.  ). b  - 

subregionC . last. ). b)  + 
sqr( subregionC  .  new.  ).  c  ~ 

subregionC . last. ). c)  ); 
if  temp  <=  tabc_dif ference  then  join  :=  true 
else  join  :=  false; 

end;     C*  sub  procedure  compareabc  *) 

procedure  checksdCfirst , second: coefrec;  var  join: boolean); 

var 
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ii,jj  :  integer; 
begin 

for  ii  : =  1  to  3  do 
for  j j  :=  1  to  4  do 

f  irst.  constarr( .  ii,  j  j.  )  :  = 

f irst. constarr(. ii, j j. )  + 
second. constarr(. ii,jj. ); 
first.ee  :=  first.ee  +  second.ee  ; 
first,  ex  :=  first,  ex  +  second,  ex  ; 
f  irst.  ey  :=  f  irst.  ey  +  second,  ey  ; 
with  first  do 
begin 

gauss(constarr,a,b,c); 

sd       :=  sqrt(variance( a,b,c,ee,ex,ey jconstarr) ); 
if  sd  >  max_sd  then  join  :=  false  else  join  :=  true  ; 
end; 
end;    (*  sub  procedure  check_sd_before_combine  *) 

begin       (*   procedure  combine  *) 
tabc_dif f erence  :  =  increment; 

while  tabc_dif ference  <=  abc_dif ference  do  begin 
joined  : =  false; 
for  i  :=  2  to  (row-2)  do 
begin 

last  :=  pt(.  i,2.  ).  gp  ; 
for  j  :=  2  to  (col-2)  do 
begin 

new  :  =  pt(.  i,j  +  l.  ).  gp; 

if  ((i=row-2)  and  (j=col-2))  and  (not  joined)  then 

subregion( .  new.  ).  active  :=  true; 
if  (newolast)  and 

(subregion( . new. ). maxmag  <=  upper_bound)  and 
( subregion( . last.  ) . maxraag  <=  upper_bound)  then 
begin 

compareabc( joinabc); 
if  joinabc 
then  begin 

checksd(  subregion( . last. ) , 

subregion( . new. ) , joinsd  ); 
if  joinsd  then  begin 

do_combine( pt , 

subregion, last , new); 
joined  : =  true; 
end 
else  begin 

subregion(.  last.  ).  active  :=  true; 
joined  :  =  false  ; 
last  : =  new  ; 
end; 
end 
else  begin 

subregionC .  last.  ).  active  :=  true; 
joined  :=  false  ; 
last  : =  new  ; 
end; 
end 
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else  begin 

subregion( .  last.  )•  active  :=  true; 
joined  :=  false  ; 
last  : =  new  ; 
end; 
end; 
end;     (*   row  combine  finished  *) 
for  j  : =  2  to  (row-2)  do 
begin 

last  :=  pt(.  2,  j.  ).  gp  ; 
for  i  :=  2  to  (col-2)  do 
begin 

new  :=  pt(.  i+1,  j.  ).  gp; 
if  (newolast)  and 

(subregion(.  new.  )•  niirirnsg  *-=  upper_bound)  and 
(subregion(.  last.  ).  rainmag  <=  upper_boimd)  then 
begin 

compareabc( joinabc); 
if  joinabc 
then  begin 

checksd(  subregion( . last. ) , 

subregion( . new. ) , joinsd  ); 
if  joinsd  then  begin 

do_corabine(pt , 

subregion, last ,new); 
joined  : =  true; 
end 
else  begin 

subregion( .  last.  ).  active  :=  true; 
joined  :=  false  ; 
last  :  =  new  ; 
end; 
end 
else  begin 

subregion( .  last.  ).  active  :=  true; 
joined  : =  false  ; 
last  : =  new  ; 
end; 
end 
else  begin 

subregion(.  last.  ).  active  :=  true; 
joined  :  =  false  ; 
last  : =  new  ; 
end; 
end; 
end;     (*  column  combine  finished  *) 
tabc_dif f erence  :  =  tabc_dif f erence  +  increment; 
end;   (*  while  ''O 

end;      i*   procedure  combine  *) 

procedure  cal_avg_dev(var  subregion: subregion_array; 

var  points: point_array;  var  avg: real); 
var 

sum:  real; 
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i,  j,gp:  integer; 
begin 

sum  :  =  0.0; 

for  i: =  2  to  row-1  do 

for  j:=  2  to  col-1  do 

begin 

gp  :=  points (.  i,j.  )•  gp; 

sum  :=  sura  +  abs((  subregion( .  gp.  )•  s*i  ■•■ 

subregion(  .  gp.  ). b*j  + 
subregion(.  gp.  ).  c) 
-  points(. i, j. ), el); 
end; 

avg  :=  sum  /  ( (  row-2)''f(col-2)  ); 
writeln( ' average  deviation  =  ' ,avg  :  15); 
end;  {*   sub  procedure  cal_average_deviation  *) 


procedure  writedata(points:  point_array); 
var 

r,c  :  integer; 
tempa,tempb, tempo  :  real; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 
begin 

tempa  :=  subregion( .  points( .  r  ,c.  ).  gp.  ).  a; 
tempb  :=  subregion( .  points( .  r  ,c.  ).  gp.  ).  b; 
tempo  :=  subregion( .  points( .  r , c.  ) .  gp.  ).  c; 
points( .  r  ,c.  ).  el  :=  trunc  (tempa''<'r  +  tempb''«'c  +  tempo); 
end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

writeln(doutl ,points( .  r,c.  ).  el); 
end; 

('■'■  procedure  writedata  *) 


procedure  writetempdata(points: point_array); 
var 

k, r , c,rl , r2 , cl , c2  :  integer; 

tempa, tempb, tempc  :  real; 

current  :  nodepointer; 
begin 

for  r: =  2  to  row-1  do 

for  c: =  2  to  col-1  do 

points( .  r  ,c.  ).  el  :=  1  ; 

for  k: =  1  to  counter  do 

if  subregion(.  k. ). active  then 

begin 

current  :=  subregion( . k. ). list  ; 
while  current  <>  nil  do  begin 
rl:  =current-'.  r; 
cl:  =current-'.  c; 

if  current-",  next  <>  nil  then  begin 
r2:  =current-'.  next--,  r; 
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c2:  =current-'.  next--,  c; 
end 

else  begin 
r2  :=  rl; 
c2  :=  cl; 
end; 
if  rl  >  r2  then  begin  r:  =  r2;  r2:  =  rl;  rl:=  r;  end; 


if  cl  >  c2  then  begin  c:  =  c2;  c2:  =  cl;  cl:  =  c;  end; 

tempa  :=  subregion( .  k.  )•  s; 

tempb  :=  subregion( .  k.  )•  b; 

tempc  :=  subregion( .  k.  )•  c; 

if  rl=r2  then  for  c  :=  cl  to  c2  do 

points( .  rl  ,c.  )•  el  '•-   trunc  (tempa^rl  +  tempM'c  +  tempc); 

if  cl=c2  then  for  r  :=  rl  to  r2  do 

pointsC . r ,cl. )• el  :=  trunc  (tempa"r  +  tempb^cl  +  tempc); 

current  :  =  current--,  next; 
end; 
end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

writeln(dout2,points( . r ,c. ). el); 
end; 

("  procedure  writetempdata  *'^) 


procedure  check; 
var 

num,r,c,i  :  integer; 
begin 

for  c:  =  1  to  counter  do  temparr(.c.  )  :=  false;  (-  check  -) 

num  :  =  0; 

for  i  : =  1  to  counter  do 

if  subregion(. i. ). active  then 
begin 

r:=  i  ; 

temparr(.r.  )  :=  true  ; 

num  : =  num  +1  ; 

while  r  <>  0  do 

begin 

with  subregion( . r. )  do 
temparr(.r.  )  :=  true; 
r:  =  subregion( .  r.  ).  side  ; 
end; 
end; 

for  c:  =  1  to  counter  do  temparr(.c.  )  :=  false;  (*  check  *) 

num  :  =  0; 

for  i  : =  1  to  counter  do 

if  subregion(.  i. ).  active  then 
begin 

r:=  i  ; 

temparr(.r>)  :=  true  ; 

num  : =  num  +1  ; 
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while  r  <>  0  do 
begin 

with  subregion( .  r.  )  do 
temparr(.r.  )  :=  true; 
r:  =  subregion( .  r.  )•  side  ; 
end; 
end; 

writeln('num  of  gps  =  ', counter  :3,'  *  lowbound=' , low_bound: 3, 

'  upper_bound=' ,upper_bound:  3, '   '"f  tnax_sd=' ,max_sd:  3  ); 
writeln( 'unexamined  gps  are  =>'); 
c:=  0; 
for  i  : =  1  to  counter  do  if  temparr(.i. )  =  true  then  c: =  c+1 

else  writeC*'  ,i:  3, ''>'); 
writeln; 

writeln('#  of  gps  examined  =>',c:4,'  *  abc_difference=' , 
abc_dif ference); 

writeln('//  of  combined  gps  =>' ,num  :4); 
end;    (''=■  procedure  check  •") 
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Main  Program 


begin 
page; 

initialize  ; 

new( top); 

top--,  up  :  =  nil  ; 

reset(data, 'data  ell'); 

reset(datam, ' data  magi'); 

rewrite(doutl, 'data  outl');    (*  data  after  combine  process  *) 

rewrite(dout2, ' tdata  outl');   (*  data  of  boundary  *) 

for  c  :  =  1  to  col  do 
for  r  : =  row  downto  1  do 

readln(data,points(. r. )(. c. ). el); 
for  c  : =  1  to  col  do 
for  r  :  =  row  downto  1  do 

r eadln(dat am, point s( .  r.  )( . c. ). mag); 

makesub(2,(row-l) ,2,(col-l) ,top); 

printsubl(points);  (*   print  the  subregions  created  by  makesub  *) 

for  i  :  =  1  to  counter  do 
with  subregionC . i. )  do 
begin 

get const (points ,crl,crh,ccl ,ccr,ok,ee,ex,ey); 

gauss(ok,a,b,c); 

sd       :=  sqrt(variance(a,b, c,ee,ex,ey ,ok) ); 

i*   sd  for  region  i  after  division  *) 

constarr  :  =  ok  ; 
end; 
combine (points , subregion); 
get_nodes( subregion) ; 

printsub2(points);  {*   print  the  subregions  created  by  combine  ''«■) 
writedata(points)  ; 

i*   write  modified  elevation  data  into  file  'data  out/Z'^'O 
writetempdata( points); 

i*   write  elevation  data  into  file  'tempdata  out//'  ''<■) 
cal_avg_dev( subregion, points , avg) ; 
i*   check;  *) 
end. 
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APPENDIX  B.     PASCAL  CODE  OF  STRICT  BOTTOM-UP  ALGORITHM 


(^^$s450000'O 


3. 


Author   =  Yee,  Seung  Hee. 
Date     =  3  June  1988. 

Input    =  1.  Elevation  data  acquired  from  Fort.  Hunter  Liggett 
data  base  (40  by  40  square). 
2.  file  of  magnitude  of  gradient  of  the  above  elevation 
data  created  by  using  the  program  in  Appendix  3. 
file  of  curvature  of  points  created  by  using  the 
program  in  Appendix  3. 

Elevation  data  of  the  model  created  by  this  program. 
Elevation  data  of  the  boundary  points,  other  points 
are  given  the  value  one  so  that  when  we  display 
this  model  in  three-dimensional  pictures,  we  can 
only  see  the  boundaries  of  the  patches. 
Computer  =  Waterloo  Pascal  on  IBM  370/3033AP, 

Naval  Postgraduate  School. 
Description  = 

This  program  is  the  Pascal  implementation  of  the 
Joint  top-down  and  bottom-up  terrain  modeling. 


Output 


1. 
2. 


J^  J^  k*^  y .  J^  >'^  ^.  J^  JU  J^  ^^  ^^  mS^  , 


'^^^J.J^J.. 


'fVf*VfV.- 


i<-hi(  • 


PROGRAM  Bottom_up( input , output ) ; 
const 

low_bound  =  0. 1;    ("  magnitude  <  10  %  then  level  slope 

("  magnitude  >  50  ?o  then  unsafe  slope 


up_bound  =  0.  6;     (" 
magnitude_ratio  =  0.  2;    (" 


else  safe  slope 
ratio  of  magnitude  : 

dif  of  two  mags  /  old  mag 


abc_comparison_const  =20; 
max_group_num  =  400; 
row  =  40; 
col  =  40; 
type 

elrange  =  1.  .  10000; 

mag_type  =  ( level, safe, unsafe); 

point   =  record 

r  ,c:  integer; 
end; 
nodepointer  =  --node; 
node  =  record 

r ,c  : integer; 
next  :  nodepointer; 
end; 
point_rec  =  record 

el  :  elrange  ; 

mag  :  real; 

magtype:  mag_type; 

cur:  char; 

gp  :  0. . max_group_num; 

ptnext  :  point; 


Vc) 
") 
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var 


end; 
point_array  =  array  ( .  1.  .  row,  1.  .  col.  )  of  point_rec; 
polynomial_array  =  array  ( .  1.  .  3, 1.  .  4,  )  of  real; 
coef_record  =  record 

start  :  point; 

upper_left  :  point; 

a,b,c  :  real  ; 

sd   :  real  ; 

ee,ex,ey  :  real; 

side  :  0. . max_group_num; 

poly_array  :  polynotnial_array; 

active  :  boolean; 

list   :  nodepointer  ; 

numpts  :  integer; 

avgmag:  real; 
end; 
coef_array  =  array  (. 1.  . max_group_num.  )  of  coef_record  ; 

points  :  point_array  ; 
subregion   :  coef_array  ; 
ok     :  polynomial_array; 

data,datara,datac,doutl,dout2,dout3,dout4  :  text; 
counter ,r,c, i,num  :  integer; 
tempmag  :  real; 


function  variance  (a,b,c,ee,ex,ey: real;  ok: polynomial_array):  real; 

begin 

variance  :=  (ee-(  Z^'^a^^ex) -(  2*b'''-ey) -(  2''''c''^ok( .  3,4.  ) ) 
+( a'''a*ok( .  1 , 1.  ) )+( 2->a^^b''^ok(  .2,1.)) 
+( 2'-b'Vc'Vok( .  2 , 3.  ) )+( 2*a'Vc*ok(  .1,3.)) 
+(b'>b''^ok( .  2 , 2.  ) )+( c*c'^ok(  .3,3.))) 
/(ok(.3,3.  )-l); 


end; 

proce 

dure  initialize; 

var 

r, 

c  : integer; 

begin 

counter  : = 

0  ; 

for  c  :  =  1 

to  co] 

do 

for  r  :=  1 

to  ] 

row  do 

begin 

points(. 

r,c. 

)• 

gP  :=  C 

; 

points(. 

r,c. 

). 

ptnext. 

r  :  = 

0  ; 

points(. 

r,c. 

)• 

ptnext. 

c  :  = 

0  ; 

end; 

for  r  : =  1 

to  max 

:_group_ 

.num  do 

begin 

subreg. 

Lon( . 

r. 

).a  :  = 

-999; 

subreg. 

Lon( . 

r. 

).b  :  = 

-999; 

subreg. 

Lon( , 

r. 

).c  :  = 

-999; 

subreg. 

Lon( 

r. 

). active  : = 

true; 

subreg. 

Lon( 

r. 

).  side 

:=  0; 

subreg. 

Lon( 

r. 

).sd  :  = 

=  -999; 

end; 


end; 


(*   procedure  initialize  *) 
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procedure   getconst(  r,c: integer; 

var  arr:  polynomial_array;  var  ee,ex,ey  : real); 


begin 


arr(.l,l.  ) 
arr(.l,2.  ) 
arr(. 1,3.  ) 
arr(.l,4.  ) 

arr(.2,l.  ) 
arr(. 2,2. ) 
arr(.2,3.  ) 
arr(.2,4.  ) 

arr(.  3,1.  ) 
arr(. 3,2.  ) 
arr(. 3,3.  ) 
arr(.3,4.  ) 


=arr(.  1,1.  )+(r*r); 

=arr(.l,2.  )+(c'Vr); 

=arr(.  1,3.  )+  r; 

=arr(.  1,4.  )+(points(. r,c. ). el*r); 

=arr(.2,l.  )+(r*c); 

=arr(.2,2.  )+(c*c); 

=arr(.  2,3.  )+  c; 

=arr(.  2,4.  )+(points(.  r,c.  ).  el*c); 

=arr(.  3,1.  )+r; 
=arr(.  3,2.  )+c; 
=arr(.  3,3.  )+l; 
=arr(.  3,4.  )+points(.  r,c.  ).  el; 


end; 


ee 
ex 

ey 


=  ee  +  sqr(points( .  r  ,c.  ).  el); 
=  ex  +  (points(. r,c. ). el  *  r  ); 
=  ey  +  (points( . r,c. ). el  *   c  ); 


(*  procedure  get  const  *) 
procedure  gauss(ok: polynomial_array;  var  aa,bb,cc: real); 

var 


i,j  :  integer; 

temp  :  array  (.1..4.  )  of  real; 
t  :  real; 
begin 


('''   error  checking  purpose  ") 


i  :=  1; 

while  ok(.  i,l.  )  =  0.  0  do   i 

for  j  :=  1  to  4  do 

begin 


=  i  +  1 


(*  row  rotation  *) 


tempC  j.  ) 
ok(.l,j.) 
ok(.  i,  j.  ) 

end; 

t  :=  ok(.  1,4.  ); 


=  ok(.l,j.); 
=  ok(. i,j.  ); 
=  terap(. j.  ); 


for  i  : =  2  to  3  do 

for  j  : =  4  downto  1  do 

ok(.i,j.)  :=  ok(.i,j.)*ok(.l,l.  )  +  ok(.l,j.)*(-ok(.i,l.  )); 


i  :=  2; 

while  ok(.i,2.  )  =  0.  0  do   i 

for  j  :=  1  to  4  do 

begin 


=  i  +  1 


end; 


temp(. j. ) 
ok(.2,j.) 
ok(.  i,j.  ) 


=  ok(.2,j.); 
=  ok(. i,j.  ); 
=  temp(. j. ); 


for  j  :  =  4  downto  1  do 
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ok(.3,j.) 


ok(.3,j.)''^ok(.2,2.  )  +  ok(.2,j.)'H-ok(.3,2.  )); 


if  ok(. 3,3.  )  <>  0.  0  then 
begin 

ok(.3,4.)  :=  ok(.3,4.)  /  ok(.3,3.); 

ok(.2,4.)  :=  (  ok(.2,4.)  -  ok( .  3 ,4.  )*ok( .  2,3.  ) )  /  ok(.2,2.); 

ok(.l,4.  )  :=  (  ok(.l,4.  )  -  ok( .  2,4.  )''^ok(.  1,2.  ) 

-  ok(.3,4.  )'Vok(.l,3.  ))  /  ok(.l,l.  )  ; 
aa  :=  ok(.  1,4.  );   (*  solution  *) 
bb  :=  ok(.2,4.  );  {*   solution  *) 
cc  :=  ok(.3,4.  );  (*   solution  *) 
(*writeln( 'gauss  error  =' , 

t-(   aa''^ok(.  1,1.  )+bb'^ok(.  1,2.  )+  cc*ok(.  1,3.  )))*) 
end; 
end;      (*  procedure  gauss  *) 

procedure  prefix(nura: integer); 
begin 

case  ((nura  div  100)  mod  10)  of 


end; 


0 

write( 

1 

write( 

. 

2 

write( 

: 

3 

write( 

- 

4 

write( 

+ 

5 

write( 

% 

6 

write( 

Vc 

7 

write( 

& 

8 

write( 

7 

9 

write('=') 

end; 

(,V 

procedur 

a  prefix 

■i; 

) 

procedure  priRtsub(points: point_array); 
var 

r,c  :  integer; 
begin 
page; 

for  r  : =  1  to  row  do 
for  c  : =  1  to  col  do 
begin 

if  c<>  col  then  begin 

pref ix(points(.  r ,c.  ). gp); 
write( (points(. r,c.  ).  gp  mod  100): 2); 
end 
else  begin 

pref ix( points ( .  r ,c.  ).  gp); 
writeln((points(. r,c.  ). gp  mod  100): 2); 
writeln; 
end; 
end; 
end;      (*   procedure  printsub  *) 


procedure  printsub2(points: point_array); 

var 

k,r ,c,rl,r2,cl,c2  :  integer; 
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current  :  nodepointer; 
begin 
page; 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 

points(.  r,c.  ).  gp  :=  0  ; 
for  k: =  1  to  counter  do 
if  subregion( .  k.  )• active  then 
begin 

current  :=  subregion( .  k.  )•  list  ; 

while  current  <>  nil  do 

begin 

rl:  =current-',  r; 
cl:  =current-'.  c; 

if  current--,  next  <>  nil  then  begin 
r2:  =current-'.  next-*,  r; 
c2:  =current-'.  next-",  c; 
end 

else  begin 
r2  :=  rl; 
c2  :=  cl; 
end; 

if  rl  >  r2  then  begin  r:  =  r2;  r2:=  rl;  rl:=  r;  end; 
if  cl  >  c2  then  begin  c:  =  c2;  c2:  =  cl;  cl:  =  c;  end; 

if  rl=r2  then  for  c  :=  cl  to  c2  do 

points(.  rl,c.  )•  gp  '•-   k; 

if  cl=c2  then  for  r  :=  rl  to  r2  do 

points(.  r ,cl.  )•  gp  :=  k; 
current  :  =  current-',  next; 
end; 


end; 

for  r 

:=  1  to 

row  do 

for  c 

:=  1  to 

col   do 

begin 

if 

c<>  col 

then  b 

pref ix(points(.  r,c.  ).  gp); 
write((points(. r,c. ). gp  mod  100): 2); 
end 
else  begin 

pref ix(points( . r,c. ). gp); 
writeln( (points( .  r ,c. ). gp  mod  100): 2); 
writeln; 
end; 
end; 
end; 

(*  procedure  printsub2  *) 

procedure  printmagcur; 
var 

r,c  :  integer; 
begin 


(*  page'^) 

for  r 

: =  1  to  row 

do 

for  c 

:=  1  to  col 

do 

begin 
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if  c<>  col  then  begin 

case  points( .  r  ,c.  )•  n^sgtype  of 

level  : write( ' l' ,points(.  r ,c.  )• cur , '  '); 
safe   : write( ' s ' ,points( .  r ,c.  )• cur , '  '); 
unsafe   :write('u  ,points( . r ,c. )• cur , '  '); 
end;  (*  case  ''0 
end 
else  begin 

case  points(. r ,c. )•  roflgtype  of 

level  : writeln( ' 1' jpoints(.  r,c.  )•  cur , '  ') 
safe   : writeln( ' s' ,points( . r ,c.  )• cur , '  ') 
unsafe: writeln( 'u' , points ( .  r,c.  ).  cur, '  ' ) 
end; 

writeln; 
end; 
end; 
end;      i*   procedure  printmagcur  *) 

procedure  makeone_rowl(var  pointl,point2  :  point); 
begin 

points( .  pointl.  r  ,pointl.  c.  )•  ptnext.  r  :=  point2.  r; 
points( .  pointl.  r,pointl.  c.  )•  ptnext.  c  :=  point2.  c; 
points  ( .  point 2.  r, point 2.  c.  )•  gP  '■- 

points ( .  pointl.  r , point  1.  c.  ).  gp; 
with  subregion( . counter.  )  do 
nurapts  : =  numpts  +1; 


end;      (*  sub  procedure  make  one_l  ''<■) 


procedure  makeone_row2(var  pointl, point2  :  point); 
begin 

points( .  pointl.  r, pointl.  c.  )•  ptnext.  r  :=  point2.  r; 
points( .  pointl.  r, pointl.  c.  ).  ptnext.  c  :=  point2.  c; 
point s(  .  point 2.  r, point 2.  c.  ).  gp  :  = 

point s( .  pointl.  r, pointl.  c.  ).  gp; 
with  subregion( .  counter.  )  do 
begin 

numpts  :  =  numpts  +1; 

avgmag  :=  avgmag''''(numpts-l)/numpts  + 

points(. point2.  r,point2.  c.  ). mag/numpts; 
end; 

end;      (-  sub  procedure  make  one_2  *) 


procedure  rowlink; 
var 

i,j  :  integer; 

point l,point2  :  point; 
begin       ("  sub  procedure  rowlink  *) 

counter  :  =  0; 

for  i  : =  2  to  (row-1)  do 

begin 
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counter  : =  counter  +  1  ; 
points( .  i  ,2.  )•  gP  '•-   counter; 
subregion( .  counter.  )•  start,  r  :=  i; 
subregion( .  counter.  ).  start,  c  :=  2; 
subregionC .  counter.  ).upper_left  :  = 

subregion( .  counter.  ).  start  ; 
subregion( . counter. ). avgmag   :=  points( . i,2. ). mag; 
subregionC .  counter.  ).  numpts  '•  —   ^   > 
for  j  : =  2  to  (col-2)  do 
begin 

pointl.  r  :  =  i; 

pointl.  c  :  =  j; 

point2.  r  :  =  i; 

point2.  c  :  =  j+1; 

if  ((points( . i, j+1. ).  magtype  <>  safe)  and 

( points (. i,j. ). magtype  =  points(. i, j+1. ). magtype))  then 
makeone_rowl( pointl ,point2) 

('''  (level  level)  or  (unsafe  unsafe)  *) 
else  if  ( (points( . i,j. ). magtype  =  safe)  and 

(points( . i, j+1. ). magtype  =  safe))  and 
(points(. i, j. ). cur  =  points(. i,j+l. ). cur)  then 
if  (abs(points( .  i,  j  +  1.  ).  mag- 

subregion( .  counter.  ) .  avgmag) / 
subregion(.  counter.  ). avgmag)  < 
magnitude_ratio  then 
makeone_row2( pointl ,point2) 
else  begin 

counter  : =  counter  +1  ; 
points(.  i,  j  +  1.  ).  gp  :=  counter; 
subregion( . counter. ).  start  :=  point2; 
subregion( . counter. ). upper_left  :=  point2  ; 
subregion(. counter. ). avgmag  := 

points ( .  i,  j  +  1.  ).  mag; 
subregion( . counter.  ). numpts  :=  1  ; 
pointl  :=  point2; 
end 

else  begin 

counter  :  =  counter  +1  ; 
points( .  i,  j  +  1.  ).  gp  :=  counter; 
subregion( . counter. ). start  :=  point2; 
subregion( . counter. ). upper_left  :=  point2  ; 
subregion( . counter. ). avgmag  := 

points(.  i,j+l.  ).mag; 
subregionC. counter. ). numpts  :=  1  ; 
pointl  :=  point2; 
end; 
end; 
end; 
end;       (*  sub  procedure  row  link  *) 


procedure  makeone_col( last , new  :  point); 
var 


ii,jj  :  integer; 
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tail, current   :  point; 

begin 

subregion(.  points( .  new.  r, new.  c.  )•  gP-  )•  Active  :=  false  ; 
with  subregion(. points( . last. r, last.  c.  ). gp.  )  do 
begin 

numpts: =numpts  + 

subregion(. points(.  new. r,new.  c.  ).  gp.  ). numpts; 
avgmag:  =avgmag''f(  numpts - 

subregion(. points(. new.  r,new.  c.  ).  gp.  ).  numpts)/ 
numpts  + 

subregion( . points ( . new.  r,new.  c.  ).  gp.  ).  avgmag* 
subregion( . points( .  new.  r ,new.  c.  ).  gp.  ). numpts/numpts; 
end;       {*   reassigning  number  of  points  and  average  magnitude  ■*) 

if  (subregion( . points( . new. r,new.  c.  ).  gp. ).  upper_left.  r  + 
subregion(. points(. new. r,new. c.  ).  gp.  ).  upper_left. c)  < 
(subregion( . point s( . last. r, last. c. ). gp. ). upper_left. r  + 
subregion( . points ( . last. r, last. c.  ). gp. ). upper_left.  c)  then 

subregion( . point s( . last. r , last. c. ). gp. ). upper_left  : = 
subregion( . points( . new. r ,new.  c.  ). gp. ). upper_left  ; 

current  :=  points(. last. r, last. c. ). ptnext  ; 

if  current,  r  =  0  then  (*   or  current,  c  =  0  ">'') 

begin 

point s( .  last.  r,last.  c. ). ptnext   :  = 

subregion(. points (. new. r,new.  c.  ). gp.  ).  start  ; 
current  :=  subregion( . points( . new.  r, new.  c.  ).  gp.  ). start  ; 
end 
else  begin 

while  current. r  <>  0  do 
begin 

tail  : =  current; 

current  :=  points( . current. r, current. c, ). ptnext; 
end;  {">''     traverse  ptnext  of  lastpt  link 

to  the  end  of  link  *) 
points( .  tail,  r, tail.  c.  ).  ptnext:  =   (''''  now  link  the  tail 

points  to  the  head  of'^) 
subregion( .  points ( .  new.  r ,new.  c.  ). gp. ). start  ; 
i*   newpt  link  *) 
current  :=  (*   into  existing  link   ''''■) 

subregion( . points (,  new.  r ,new.  c.  ). gp.  ). start  ; 
end; 
while  current. r  <>  0  do  (*  reassigning  gp  #  to  new  pts  *) 
begin 

points( .  current,  r,  current,  c.  ).  gp  :=  points( .  last,  r,  last.  c.  ).  gp; 
current  :=  points( . current. r, current. c. ). ptnext; 
end;   (*  while  -0 
end;      (''=■  procedure  makeone_col  *) 


procedure  collink; 
var 

i,j  :  integer; 

point ljpoint2  :  point; 
begin      ("  sub  procedure  col  link  *) 
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for  j  :=  2  to  (col-1)  do 
begin 

for  i  : =  2  to  (row-2)  do 
begin 

pointl.  r  :  =  i; 
pointl.  c  :  =  j; 
point2.  r  :  =  i+1; 
point2.  c  :  =  j; 

if  (points(. i+1, j. )• gp  <>   points( . i, j. )• gp)  then 
if  (( point s(. i+1, j, ). raagtype  <>  safe)  and 

( points (. i,j. )• magtype  =  points(. i+1, j. )•  magtype) )  then 
makeone_col(pointl ,point2) 
i*   ( level, level)  or  (unsafe, unsafe)  *) 
else  if  ( (points (. i,j. )• magtype  =  safe)  and 

( point s(.  i+1,  j.  )•  nifigtype  =  safe))  then 
if  (points( . i, j. )• cur  =  points(. i+1, j. )• cur)  then 
if  abs(subregion( . points(.  i+1, j.  ).  gp.  ).  avgmag- 
subregion( . point s( . i, j.  ).  gp.  ). avgmag)/ 
subregion( . points ( . i , j .  ) .  gp.  ) . avgmag  < 
magnitude_ratio 
then 

niakeone_co  1  (  point  1 ,  point  2  ) ; 
pointl  :=  point2; 
end; 
end; 
end;  (''''  sub  procedure  col  link  *) 

procedure  do_combine( last , new  :  point); 

var 

ii,jj  :  integer; 
tail, current   :  point; 

begin 

subregion(. points( . new. r, new. c. ). gp. ). active  :=  false  ; 
with  subregion( . points(.  last. r, last.  c.  ).  gp.  )  do 
begin 

numpts: =numpts  + 

subregion( . points( . new, r ,new. c. ).  gp.  ).  numpts; 
end;       (■''  reassigning  number  of  points  *) 

if  (subregion(.  points( .  new.  r  ,new.  c,  ).  gp.  ).  upper_left.  r  + 
subregion( . points( . new. r ,new.  c. ).  gp.  ).  upper_left.  c)  < 
( subregion( . points( .  last,  r , last.  c.  ) .  gp.  ) .  upper_lef t.  r  + 
subregion( . points ( . last,  r , last.  c.  )• gp. )■ upper_left.  c)  then 

subregion( . points ( . last. r, last. c. ). gp. ). upper_left  :  = 
subregion( . point s ( . new. r , new. c. ) . gp. ) . upper_lef t  ; 

current  :=  points( . last. r, last. c. ). ptnext  ; 

if  current. r  =  0  then         (*   or  current. c  =  0  *) 

begin 

points( . last. r, last.  c.  ).  ptnext  :  = 

subregion( .  points( .  new.  r  ,new.  c.  ).  gp.  ).  start  ; 

current  :=  subregion(. points(. new. r, new.  c.  ).  gp.  ).  start  ; 
end 
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else  begin 

while  current,  r  <>  0  do 
begin 

tail  : =  current; 

current  : =  points( . current. r, current,  c.  )•  ptnext; 
end;  ('''  traverse  ptnext  to  link  the  head  *) 

pointsC .  tail,  r,  tail.  c.  ).  ptnext:  =  (''"  of  the  ''=■) 

subregion(. points(. new. r,new. c. ). gp. ). start  ; 
(*   new  point  *) 
current  :=  i*   into  existing  link  *) 

subregion(. points(.  new.  r,new.  c.  ).  gp.  ).  start  ; 
end; 
while  current. r  <>  0  do  (*  reassigning  gp  #  to  new  pts  *) 
begin 

points( .  current,  r, current,  c.  ).  gp  :=  points( .  last,  r,  last.  c.  ).  gp; 
current  :=  points( . current. r , current. c.  ).  ptnext; 
end;   (''f  while  *) 


poly_array(. ii,jj. )  := 


for  ii  : =  1  to  3  do 
for  j j  :  =  1  to  4  do 
subregion( .  point s( .  last,  r , last,  c.  ).  gp 

subregion( .  points( . last. r , last. c 

subregion( . points( . new. r ,new. c. ) 
subregion( . point s( .  last,  r , last.  c.  ).  gp 

subregion( . point s( . last. r , last. c 

subregion( . points( . new. r ,new. c. ) 
subregion( . point s( .  last. r , last. c. ). gp 

subregion( . points( . last. r , last. c 

subregion( . points( . new. r ,new. c. ) 
subregion( . points( .  last,  r , last.  c.  ).  gp 

subregion( . points ( . last. r , last. c 

subregion( .  points( . new. r ,new. c.  ) 
with  subregion( . points( .  last. r , last,  c 
begin 

sd       :=  sqrt( variance( a,b,Cjee,ex,ey ,poly_array) ); 
end; 
end;      ('='=■  procedure  do_combine  *) 


gp 
gp-  ) 
ee 

gp 
gp-  ) 
ex 

gp 
gp-  ) 
ey 

gp 
gp-  ) 
)-gP 


). poly_array(. ii, jj.  )  + 
poly_array(.  ii,jj.  ); 

).  ee  + 
ee  ; 

).  ex  + 
ex  ; 

).ey  + 
ey  ; 
)  do 


procedure  combine(var  pt:  point_array; 

var  subregion: coef_array); 
var 

i,j  :  integer; 

last, new  :  point; 

joinabc,joinsd, joined  :  boolean; 

procedure  compareabc(var  join:  boolean); 
var 

temp  :  real  ; 
begin 

temp  :=  sqrt(  sqr(subregion( . points( . new.  r ,new. c. ).  gp.  ).  a 

subregion(. points (. last.  r,last.  c.  ).  gp.  ).  a) 

sqr(subregion( . points( . new. r ,new. c. ). gp. ).  b 

subregion(.  points (.  last.  r,last.  c.  ).  gp.  ).  b) 

sqr(subregion( . points( . new. r ,new. c. ). gp.  ).  c 
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subregion( . points ( . last. r , last. c.  )•  gp-  )•  c)  ); 
if  temp  <=  abc_comparison_const  then  join  :=  true 
else  join  :=  false; 

end;      {*   sub  procedure  compareabc  *) 


begin       {'>''   procedure  combine  *) 
joined  :=  false; 
for  i  :=  2  to  (row-2)  do 
begin 

last,  r  :  =  i; 

last,  c  :  =  2; 

for  j  :=  2  to  (col-2)  do 

begin 

new.  r  :  =  i; 

new.  c  :  =  j+1; 

if  ((i=row-2)  and  (j=col-2))  and  (not  joined)  then 

subregion( . points ( . new. r ,new.  c.  ).  gp.  ).  active  :  = 

true; 
if  (points(. new. r ,new. c. ). gp  <>  points(. last. r, last. c. ). gp) 
then  begin 

compareabc( joinabc); 
if  joinabc 
then  begin 

do_combine( last ,new); 
joined  :  =  true; 
end 
else  begin 

subregion(. points(. last. r,last.  c.  ).  gp.  ).  active  :  = 

true; 
joined  :=  false  ; 
last  : =  new  ; 
end;     i*   joinabc  *) 
end 
else  begin 

subregion( .  point s( . last. r , last. c. ). gp. ). active  :  = 

true; 
joined  :=  false  ; 
last  : =  new  ; 
end;     (->''    last  <>  new  ''<■) 
end;  (*  end  j  ''') 
end;     (*  end  i  •>'^)(*   row  combine  finished  *) 
for  j  :=  2  to  (row-2)  do 
begin 

last,  r  :=  2; 

last. c  : =  j; 

for  i  :=  2  to  (col-2)  do 

begin 

new.  r  :  =  i+1; 

new.  c  :  =  j  ; 

if  (points(.  new.  r,new.  c.  ).  gp  <>  points(. last. r, last.  c. ).  gp) 

then  begin 

compareabc( joinabc); 
if  joinabc 
then  begin 
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do_combine( last ,new); 
joined  : =  true; 
end 
else  begin 

subregion( . points ( . last. r , last. c. ). gp. ). active 

true; 
joined  :=  false  ; 
last  : =  new  ; 
end;     {*   join  abc  *) 
end 
else  begin 

subregion(. points (. last. r, last. c. ).  gp. ). active  :  = 

true; 
joined  :=  false  ; 
last  : =  new  ; 
end;     (*  last  <>  new  *) 
end; 
end;     (''='  column  combine  finished  *) 

end;      {*   procedure  combine  *) 


procedure  get_boundaries(var  subregion:  coef_array); 
type 

direction  =  (  up,down,right , left  ); 
var 

i,j,k,r,c      :  integer; 

current , head, p  :  nodepointer; 

temp  :  direction  ; 

procedure  get_next(var  p  :  nodepointer; 

var  temp  :  direction); 
var 

i,j  :  integer; 

done  :  boolean  ; 


procedure  do_dir(var  done:  boolean  ;  dir: direction  ); 
begin 

new(p); 
P-.  r  :=  i  ; 
P"-  c  :=  j  ; 
done  : =  true; 
case  dir  of 

up    :   temp  : =  up; 
down   :   temp  :  =  down; 
right  :   temp  : =  right; 
left   :   temp  : =  left  ; 
end;  i'^"   case  ''•') 
end;  (*   sub_sub_  procedure  do_dir  *) 


begin  (*  sub  procedure  get  next  *) 
done  ;=  false  ; 
i  :=  P-.  r; 
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J  ■=  V-  c; 
case  temp  of 

right  :  begin 

while  not  done  do  begin 

if  points(. i-1, j.  )•  gp  =  k 

then  do_dir(done,up)    (*  go  up  and  done  *) 
else  if(  (subregion(.  k.  )•  list--,  r  =  i) 

and  (subregion(.  k.  )•  list--,  c  =  j))  then 
begin  (*  return  to  org  point  *) 

new(p); 
P-.  r  :  =  i  ; 
p-'-c  :=  j  ; 
done  : =  true; 
end 
else  if  (points(.  i, j  +  1. )•  gp  =  k) 

then  j:=  j+1  ('"=•  go  right  *) 

else  if  ( points (.  i+l,j.  )•  gp=k  ) 

then  do_dir( done, down)  (*  go  down  and  done  *) 
else  do_dir(done,  left);     (*  go  back  and  done  ''0 
end;  {*   while  ") 
end; 
left   :  begin 

while  not  done  do  begin 

if  points( . i+1 , j.  )•  gp  -   k 

then  do_dir( done, down)  (''f  go  down  and  done'^) 
else  if  (points( .  i,  j -1.  )•  gp  =  k) 

then  j  :=  j-1  i*  go   left  •>'=■) 

else  if  (points( . i-1, j. )• gp  =  k) 

then  do_dir(done  ,up)    {*   go  up  and  done  '"f) 
else  do_dir(done,  right);  (*   go  back  and  done  ■>'") 
end;  (*  while  *) 
end; 
up    :  begin 

i  :=  i-1; 

while  not  done  do  begin 

if  points(.  i,  j-1.  ).  gp  =  k 

then  do_dir(done, left)  (*  go  left  and  done  *) 
else  if  (points( .  i-1, j. )• gp  "  k) 

then  i  :=  i-1  (''^  go  up  '*') 

else  if  (points( . i, j+1. )• gp  =  k  ) 

then  do_dir(done,right)('^  go  right  and  done*) 
else  do_dir( done, down);    (*  go  down  and  done  *) 
end;  {->''   while  ■>'') 
end; 
down   :  begin 

i  :=  i+1; 

while  not  done  do  begin 

if  points ( .i,j+l.).gp  =  k 

then  do_dir( done, right )(*  go  right  and  done*) 
else  if  points( . i+1, j. ). gp  =  k 

then  i  :=  i+1  (*  go  down  *) 

else  if  points( . i, j-1. ). gp  =  k 

then  do_dir(done, left)(*  go  left  and  done*) 
else  do_dir(done,up);      (*  go  up  and  done  *) 
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end;  {*   while  *) 
end; 
end;  {*   case  ''') 
end;  {*   sub  procedure  get_next  *) 


begin   (''f  procedure  get  boundaries  *) 
for  k  : =  1  to  counter  do 
if  subregion(. k. )• active  then 
begin 

new(p); 

p-'.  r  :  =subregion(  .  k.  )•  upper_left.  r  ; 
p-".  c  :  =subregion( .  k.  )•  upper_left.  c  ; 
while  (  points( .  p-".  r-ljp--.  c.  )•  gP  ~   k.)  do 

p--.  r  :=  p-".  r-1;    {*   move  to  the  upper  most  row  of  the  gp  *) 
while  (  points(.  p-".  rjp--.  c-1.  )•  gp  =  k)  do 

p-'.  c  :=  p-'.c-l;    (*  move  to  the  left  most  col  of  the  gp  *) 

subregion( .  k.  ).  list  :=  p  ; 
if  points( .  p--.  r  jp--.  c+1.  )•  gp  —   k  then 

begin  {*   starting  toward  right  groups  *) 

current  : =  p  ; 
temp  : =  right  ; 
repeat 

get_next ( p , t emp) ; 
current-",  next  :=  p  ; 
current  :  =  p  ; 
until  (subregion( .  k.  ).  list-",  r  =  p-*.  r) 

and  (subregion( .  k.  )•  list-",  c  =  p--.  c); 
if  (temp=down)  and  (points( .  p--.  r+1  .p--.  c.  )•  gP=k)  then 
begin  ('''■  right  and  down  points  gp  *) 

current  :=  p  ;       (*   hinged  on  start  point    ">'') 
temp  : =  down  ; 
repeat 

get_next(p, temp); 
current--,  next  :=  p  ; 
current  : =  p  ; 
until  (subregion(.  k.  ).  list--,  r  =  p--.  r) 

and  (subregion(.  k.  ).  list-",  c  =  p-'.c); 
p-'.  next  :  =  nil; 
end 
else 

p--.  next  :=  nil;  (*  right  but  not  down  gp  *) 

end 

else  if  points(.  p--.  r+1  ,p-'.  c.  ).  gp  =  k  then 
begin  ('"'  starting  toward  down  group  *) 

current  :  =  p  ; 
temp  : =  down  ; 
repeat 

get_next ( p , temp ) ; 
current--,  next  :=  p  ; 
current  : =  p  ; 
until  (subregion( .  k.  ).  list-',  r  =  p--.  r) 

and  (subregion(.  k.  ).  list-',  c  =  p-'.c); 
p--.  next  :  =  nil; 
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end 
else 

p--.  next  :=  nil;  (''=•  one  point  gp  *) 
end; 


the  following  segment  of  code  may  be  used  for  checking  purpose 

for  k  : =  1  to  counter  do 

if  subregion( . k.  )• active  then 

begin 

write(  '  Boiindary  of  group  ',k:3j'is  =>'); 
current  :=  subregion(  .  k.  )•  list; 
while  current  <>  nil  do 
begin 

write( current-",  r:  2,  '  ,  *  , current-*,  c:  2,  '=>'  ); 
current  :  =  current-",  next; 
end; 

writeln; 
end; 
■  — *  ■ 

end;   ('''■  procedure  get_boundaries  *) 

procedure  calc_coeffs; 
var 

k, i, j ,r ,c,gpnura,ptnum  :  integer; 
t,tl   :  point; 

rowchangedjColchanged:  boolean; 
begin 

gpnum  :  =  0; 
ptnum  :  =  0; 
for  k  : =  1  to  counter  do 

if  subregion( . k.  ). active  then 
begin 

rowchanged  :=  false; 
colchanged  :  =  false; 
for  i  : =  1  to  3  do 
for  j  :=  1  to  4  do 

subregion(. k. ). poly_array(.  i, j.  )  :=  0; 
subregion( .  k.  ).  ee  :=  0; 
subregion( .  k.  ).  ex  :=  0; 
subregion( .  k.  ).  ey  :=  0; 
ptnum  :=  ptnum  +subregion( . k.  ) .  nurapts  ; 

(■'<■  num  of  linked  pt  *) 
gpnum  :  =  gpnum  +1  ;  (*  num  of  active  gp  *) 

t  :=  subregion(.  k.  ).  start; 
while  t.  r  <>  0  do 
begin 

tl  :=  t; 

with  subregion( .  k.  )  do 

getconst(  t.r,  t. c,poly_array,ee,ex,ey); 
t:  =  points( .  t.  r,  t.  c.  ),  ptnext  ; 
if  (t.  r<>0)  and  (t.  r<>tl.r)  then 

rowchanged  :  =  true; 
if  (t.coO)  and  (t.  c<>tl.c)  then 
colchanged  :=  true; 
end; 
with  subregion( .  k.  )  do 
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begin 

if  rowchanged  and  colchanged  then 

gauss(poly_array ,a,b,c)         (*   normal  points  *) 
else  if  not( rowchanged)  and  not( colchanged)  then 
begin 

a  :=  0; 
b  :=  0; 

c  :=  poly_array(.  3,4.  ); 
end  i*   one  point  gp  *) 

else  if  not( rowchanged)  and  colchanged  then 
begin  (*  horizontal  line  group  *) 

a  :=  0; 

for  j:=  4  downto  2  do 
poly_array(.  3,  j.  )  :  = 

poly_array( . 3, j.  )  *  poly_array(.  2,2.  )  + 
poly_array(. 2,j. )  *   ( -poly_array( . 3,2.  ) ); 
poly_array(.  3,4.  )  :  = 

poly_array( .  3,4.  )  / 
poly_array(.  3,3.  ); 
poly_array(. 2,4. )  :=  (poly_array( .  2,4.  ) 

poly_array(.  3,4.  )  ''' 
poly_array(.2,3.  ))   / 
poly_array(. 2,2.  ); 
b  :=  poly_array(.  2,4.  ); 
c  :=  poly_array(.  3,4.  ); 
end 

else  if  rowchanged  and  not( colchanged)  then 
begin  {'''   vertical  line  group  *) 

b  :=  0; 

for  j:=  4  downto  1  do 
poly_array(.  3,j.  )  :  = 

poly_array(.  3, j.  )* 
poly_array( .  1 , 1. )  + 
poly_array(. 1, j.  )* 
(-poly_array(.  3,1.  )); 
poly_array( .  3,4.  )  :  = 

poly_array( .  3,4.  )    / 
poly_array(.  3,3.  ); 
poly_array( .  1 ,4.  )    :=   (poly_array( . 1 ,4.  ) - 
poly_array(.  3,4.  )* 
poly_array(.  1,3. ))/ 
poly_array(.  1,1.  ); 
a   :=  poly_array( .  1 ,4.  ); 
c   :=  poly_array(.  3,4.  ); 
end; 
(*  writelnC    a,b,  c,  =  ' ,a:  10 , '    ',b:10,'    '  ,c:  10)   *) 

end;  ("        with  subregion(.  k.  )  ^') 

end;  {*   if  subregion( . k. ). active  *) 

write ln( ' active  pts=' ,ptnum: 5, '  *   active  gpnum  =' ,gpnum: 3); 
end;  (''^  procedure  calc_coeffs  *) 


procedure  write_planar_data(points: point_array); 
var 

r,c  :  integer; 


terapa, tempb, tempc  :  real; 
begin 

for  r:  =  2  to  row-1  do 
for  c: =  2  to  col-1  do 
begin 

tempa  :  =  subregion( .  points(.  r  ,c.  )•  gP- )•  ^5 

tempb  :=  subregion( .  points( ,  r ,  c.  )•  gP- )  •  b; 

tempc  :=  subregion( .  points( .  r  ,c.  )•  gP- )•  c; 

points(.  r,c.  )•  el  :=  trunc  (tempa'^r  +  terapb^c  +  tempc); 
end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

writeln(doutl,points( . r,c. ). el); 
end; 

i*   procedure  write_planar_data  *) 


procedure  write_boundary_data( points: point_array); 
var 

k,r ,c,rl ,r2,cl ,c2  :  integer; 
tempa, tempb, tempc  :  real; 
current  :  nodepointer; 
begin 

for  r: =  2  to  row-1  do 
for  c: =  2  to  col-1  do 

points( .  r,c.  )•  el  :=  1  ; 
for  k: =  1  to  counter  do 
if  subregion( . k. )• active  then 
begin 

current  :=  subregion( . k.  )•  list  ; 

while  current  <>  nil  do 

begin 

rl:  =current-'.  r; 
cl:  =current-'.  c; 

if  current--,  next  <>  nil  then  begin 
r2:  =current-'.  next--,  r; 
c2:  =current-'.  next--,  c; 
end 

else  begin 
r2  :=  rl; 
c2  :=  cl; 
end; 

if  rl  >  r2  then  begin  r:  =  r2;  r2:  =  rl;  rl:=  r;  end; 
if  cl  >  c2  then  begin  c:  =  c2;  c2:=  cl;  cl:=  c;  end; 

tempa  :=  subregion( .  k.  )•  s; 

tempb  :=  subregion( .  k.  )•  b; 

tempc  :=  subregion( .  k,  )•  c; 

if  rl=r2  then  for  c  :=  cl  to  c2  do 

pointsC .  rl  ,c.  ).  el  :=  trunc  (tempa^'^rl  +  tempb*c  +  tempc); 

if  cl=c2  then  for  r  :=  rl  to  r2  do 

points( . r,cl. )• el  :=  trunc  (tempa"r  +  tempb*cl  +  tempc); 

current  :=  current--,  next; 
end; 
end; 
for  c  : =  1  to  col  do 
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for  r  : =  row  downto  1  do 

writeln(dout2 ,points( .  r,c.  ).  el); 
end;   ('"'  write_boundary_data  ''") 

(*  procedure  write_boundary_data  *) 
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Main  Program 


begin 
page; 

initialize  ; 

reset(data, 'data  ell'); 
reset(datara, ' data  magi'); 
reset(datac, ' data  curl'); 
rewrite(doutl, ' data  outl'); 
rewrite(dout2, ' tdata  outl'); 


='0 


(*  data  after  combine  process  *) 
(*  data  of  boundary  *) 


for  c 
for  r 

for  c 

for  r 
begin 


end; 
for  c 
for  r 


: =  1  to  col  do 

: =  row  downto  1  do 

readln(data,points( . r,c. ). el); 
:=  1  to  col  do 
: =  row  downto  1  do 

readln( datam , tempraag) ; 

if  tempmag  <  low_bound  then 

points( .  r  ,c.  )•  tiiagtype  :=  level 
else  if  tempraag  >  low_bound  then 

points( . r,c. )• roagtype  :=  safe 
else    points( .  r  ,c.  )•  insgtype  :=  unsafe; 
points( .  r  ,c.  )•  roag  :=  tempmag; 

: =  1  to  col  do 

: =  row  downto  1  do 

readln(datac,points( . r,c. ). cur); 


('''■printraagcur*)     (*  print  magnitude  and  curvature  *) 

rowlink;        (''=■  link  points  in  row  and  make  subregions  *) 

printsub(points);  (''"  print  subregions  created  by  rowlink  *) 

collink;        (''<■  link  points  in  row  and  column  and  make  subregions  *) 

calc_coeffs;      (*  calculate  plane  expressions  for  the 

subregions  created  ''") 
combine(points , subregion) ; 

("  combine  adjacent  regions  if  they 
have  similar  plane  expressions  *) 
calc_coeffs;  {*   calculate  plane  expressions  for  the 

final  subregions  *) 
get_boundaries( subregion); 

{*   get  boundaries  of  subregions  in  linked  list  *) 
printsub2(points);     (*   print  subregions  created  by  combine  *) 
write_planar_data( points);   (*  write  out  the  elevation  of  the  model  *) 
write_boundary_data(points); (*  write  out  the  elevation  of 

the  boundaries  of  the  subregions  *) 
end. 
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APPENDIX  C.     PASCAL  CODE  FOR  SMOOTHING  ELEVATION  DATA 

AND  CALCULATING  GRADIENT 

(*Ss300000-) 

program  calc_s lope_curv_smooth( input , output ) ; 

Author   =  Yee,  Seung  Hee. 

Date     =  3  June  1988. 

Input    =  1.  Elevation  data  acquired  from  Fort.  Hunter  Liggett 

data  base  (40  by  40  square). 
Output   =  1.  File  of  magnitude  of  gradient  of  the  points. 

2.  File  of  curvature  of  the  points  (optional). 

3.  File  of  smoothed  elevation  data  (optional). 
Computer  =  Waterloo  Pascal  on  IBM  370/3033AP, 

Naval  Postgraduate  School. 
Description  = 

This  program  contains  two  major  procedures: 
calc_slope_curvature  and  smoothing. 

The  procedure  calc_slope_curvature  calculates  the  slope 
and  curvature  of  each  points,  the  procedure  smoothing  processes  the 
elevation  data  to  get  smoothed  elevation  data. 

Original  elevation  data  is  represented  in  feet  and  there  are  40-40 
data  cells  in  one  window.  The  distance  between  each  cell  isl2.5 
meters  ,  so  it  needs  to  be  converted  to  get  a  metric  system. 

const 

distance.unit  =  1;    (^  3.  28084-12.  5  =  41.0105  -) 

row  =  3; 

col  =  3; 

curv_thresh  =  0.  05;  (•"  threshold  for  eigenvalues  of  biquadratic  ") 

upbound  =  0.  6;   (-  slope  >  60  %  then  unsafe  slope  "') 

lobound  =  0.  1;    (*  slope  <  10  %  then  level  slope  ''O 
type 

elrange  =  0. . 10000; 
pointrec  =  record 

el  :  elrange  ; 
slope:  real; 
cur:  char; 
end; 
point_array=array  (. 1. . row,l. . col. )  of  pointrec  ; 
var 

points jSpoints  :  point_array  ; 

data, outs  lope, outcur,smoutel,smouts lope, smoutcur   :  text; 

counter ,r,c,i,num:  integer; 

k2,k3,k4,k5 ,k6, slope, laml, lam2: real; 

procedure  init(var  pt:  point_array); 
var 
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r,c,i,j:  integer; 
begin 

for  r  : =  1  to  row  do 
begin 

pt(.  r,l.  ).  slope  :=  999; 

pt(.  r, col.  )•  slope  :=  999; 

pt(.  r,  1.  ).  cur  :  =  'u' ; 

pt(.  r,col.  )•  cur  :=  'u'; 
end; 

for  c  :  =  1  to  col  do 
begin 

pt(.  l,c.  ).  slope  :=  999; 

pt(.  row, c.  )•  slope  :=  999; 

pt( .  1  ,c.  ).  cur  :  =  'u' ; 

pt( .  row,c.  )•  cur  :=  'u'; 
end; 

end; 

procedure  calc_slope_curvature( var  pt:  point_array); 


var 

r,c,i,j:  integer; 

tempi , temp2 , t ermx , terray , dtermx , dtermy , 
templam   , termxx,termyy ,termxy,term  :  real; 
begin 

for  r: =  2  to  row-1  do 

for  c: =  2  to  col-1  do 


begin 

term 


=  0.  0; 


termx  :  =  0.0; 
termy  :  =  0.0; 
termxx  :  =  0.0; 
termyy  :  =  0.0; 
termxy  :  =  0.0; 


for  i  : = 

for  j  := 

begin 
term 
termx 
termy 
termxx 
termyy 
termxy 

end; 


1  to  1  do 
1  to  1  do 


pt(.  r+i,c+j. ). el; 
+  j''^pt(.  r+i,c+j.  ).el; 
=  termy  +  i''<'pt( .  r+i, c+j .  ).  el; 


=  term  + 
=  termx  + 


el; 


=  i:ermx  -^  j^^pn, .  r-^l , c-^J.  ;.  ei; 
=  termy  +  i''<'pt( .  r+i, c+j .  ).  el; 
:=  termxx  +  j''f  j''^pt( .  r+i,  c+j.  ).  ei; 
:=  termyy  +  i''fi''fpt( .  r+i,c+j.  ).  el; 
:=  termxy  +  i^j'-'^ptC .  r+i,c+j.  ).  el; 


k2  :=  termx/(6*distance_unit); 

k3  :=  termy/ (6''fdistance_unit); 

k4  :=  (termxx/2  -  term/3)/distance_unit; 

k5  :=  termxy/ (4''''distance_unit)  ; 

k6  :=  (termyy/2  -  term/3 )/distance_unit; 

laml  :=  ( ( 2'--k4+2*k6 )  + 

sqrt(  abs(  sqr(  -2''^k4-2*k6)  -4*(  2*k4'V2*k6-sqr(k5)  )  )  )  )/2; 
lam2  :=  ( ( 2''-k4+2^--k6 )  - 

sqrt(abs(sqr(-2''^k4-2'Vk6)-4*(2*k4*2*k6-sqr(k5)))))/2; 
if  abs(laral)  <  abs(lam2)  then 
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begin 

templam  :=  laml; 

laml  :=  lam2; 

lam2  : =  templam; 
end;    ("  eigenvalues  of  bigger  slopenitude  becomes  laml  '''') 
slope  :=  sqrt(sqr(k2)  +  sqr(k3)); 

pt( .  r,c.  )•  slope  :=  slope; 

if  ( laml>curv_thresh)  and  ( lam2>curv_thresh)  then 

pt(.  r,c.  )•  cur  :=  'd'     (''^  depression   *) 
else  if  ( laml<-curv_thresh)  and  ( lam2<-curv_thresh)  then 

pt(.  r,c.  )•  cur  :=  'h'      (*   hill  or  peak  *) 
else  if  (abs( laml)<=curv_thresh)  and 
(abs( lam2)<=curv_thresh)  then 
pt( .  r ,c.  )•  cur  :=  'p'     (*    planar    *) 
else  if  ( laml<-curv_thresh)  and  ( lam2>curv_thresh)  then 

pt(.  r,c.  )•  cur  :=  's'      i*  saddle    *) 

else  if  ( laml<-curv_thresh)  and 

(abs( lara2)<=curv_thresh)  then 
pt(.  r ,c.  )•  cur  :=  'r'     {*  ridge     *) 

else  if  (abs( lam2)<=curv_thresh)  and 
( laml>curv_thresh)  then 
pt( .  r  ,c.  )•  cur  :=  'v'      ("    valley    *) 
else  if  ( lam2<-curv_thresh)  and  ( laml>curv_thresh)  then 
pt( .  r  ,c.  )•  cur  :=  'a'      (*    pass      *) 
end; 
end; 


procedure  smoothing(var  points ,spoints: point_array); 
var 

r ,c, i, j , temp,dif  :  integer; 
begin 

for  r  : =  2  to  (row-1)  do 
for  c  : =  2  to  (col-1)  do 
begin 

if  (points( . r,c.  )•  slope  >  lobound) 
then  begin 

temp  :  =  0; 

for  i  :=  r-1  to  r+1  do 

temp  :=  points( .  i,c.  )•  el  +  temp; 
for  j  :=  c-1  to  c+1  do 

temp  : =  points( . r , j.  )•  el  +  temp; 
spoints(. r ,c. )• el  :=  trunc( temp/6); 
end 

else  spoints(.  r,c.  )•  el  :=  points(.  r ,c.  )•  el  ; 
end; 

for  r  :  =  1  to  row  do 
begin 

spoints(.  r,  1.  )•  el  :=  points( .  r ,  1.  )•  el; 
spoints( .  r,col.  )•  el  :=  points(.  r,col.  )•  elj 
end; 

for  c  : =  1  to  col  do 
begin 

spoints(.  l,c.  )•  el  :=  points( .  1  ,c.  )•  el; 


87 


spoints( .  roWjC.  )•  el  :=  points( .  row,c.  )•  el; 
end; 
end;   ("  procedure  smoothing  *) 


procedure  check(var  pt: point_array); 
var 

count , countb , countp , counth , countd  : integer; 
counte, count r,countv,counta  : integer; 

r,c,aajCC  :  integer; 
begin 

count  :  =0; 

countb  :  =0 

countp  : =0 

counth  : =0 

countd  :  =0 

counte  :  =0 

count r  :  =0 

countv  :  =0 

counta  :  =0 

aa  :  =  0; 

cc  :  =  0; 

for  r  : =  1  to  row  do 

for  c  : =  1  to  col  do 

if  pt( . r,c. )• slope  <  lobound  then 
aa  : =  aa  +  1 

else  if  (pt(.  r,c.  )•  slope  =  999.0)  then 
countb  : =  countb  +  1 

else  if  (pt( . r,c. ). slope  >  upbound)  then 
cc  : =  cc  +  1 

else 

case  pt( 

P   " 
d 

h 
e 

r 

V 

a 
end: 

writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
writeln( 
end; 


r  ,c.  ).  cur  of 


countp 
countd 
counth 
counte 
count r 
countv 
counta 


=  countp  +1 
=  countd  +1 
=  counth  +1 
=  counte  +1 
=  countr  +1 
=  countv  +1 
=  counta  +1 


level  points  are'  aa); 
unsafe  points  are  ,cc); 
for  safe  points: '); 
count  planar  =' , countp); 
count  hill  =' , counth); 
count  depression  =', countd); 
count  etc    =' , counte); 
count  saddle  =' , counts); 
count  ridge  =' .countr); 
count  valley  =  , countv); 
count  pass  =' , counta)*) 


( ^ ======= 

Main  Program 
^^ ^^,,^,,^^^^^^,,^^^^^^^,,^^,,^ ^3 

begin 

reset(data, 'data  ell'); 
rewrite(outslope, ' data  magi'); 
rewrite(outcur , ' data  curl  ); 
rewrite( smoutel , ' smdata  ell'); 
rewrite( smouts lope , ' smdata  raagl ' ) ; 
rewrite(smoutcur ,' smdata  curl  ); 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 

readln(data,points(.  r,c.  ).  el); 

init(points); 

init(spoints); 

calc_slope_curvature( points); 
{*   check( points)  *) 
for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 
begin 

writeln(outslope,points(.  r,c,  ).  slope); 

writeln(outcur,points(.  r,c.  ).  cur); 
end; 
for  r: =  1  to  10  do  begin 

smooth ing( points ,spoints); 

calc_s lope_curvature( spoints ) ; 

check(spoints); 

points  :  =  spoints; 
end; 

for  c  : =  1  to  col  do 
for  r  : =  row  downto  1  do 
begin 

wr it eln( smoutel , spoints (.  r,c.  ).  el); 

writeln(smoutslope,spoints( . r,c. ). slope); 

writeln(smoutcur,spoints(.  r,c.  ).  cur); 
end; 
end. 
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